LNCS 2993 



Rajeev Alur 
George J. Pappas (Eds.) 

Hybrid Systems: 
Computation 
and Control 

7th International Workshop, HSCC 2004 
Philadelphia, PA, USA, March 2004 
Proceedings 



'If Springer 




Lecture Notes in Computer Science 2993 

Edited by G. Goos, J. Hartmanis, and J. van Leeuwen 




Springer 

Berlin 
Heidelberg 
New York 
Hong Kong 
London 
Milan 
Paris 
Tokyo 




Rajeev Alur George J. Pappas (Eds.) 



Hybrid Systems: 
Computation and Control 



7th International Workshop, HSCC 2004 
Philadelphia, PA, USA, March 25-27, 2004 
Proceedings 




Springer 




Series Editors 



Gerhard Goos, Karlsruhe University, Germany 
Juris Hartmanis, Cornell University, NY, USA 
Jan van Leeuwen, Utrecht University, The Netherlands 

Volume Editors 
Rajeev Alur 

University of Pennsylvania, Department of Computer and Information Science 
3330 Walnut Street, Philadelphia, PA 19104, USA 
E-mail: alur@cis.upenn.edu 

George J. Pappas 

University of Pennsylvania, Department of Electrical and Systems Engineering 
200 South 33rd Street, Philadelphia, PA 19104, USA 
E-mail: pappasg@ee.upenn.edu 



Cataloging-in-Publication Data applied for 

A catalog record for this book is available from the Library of Congress. 

Bibliographic information published by Die Deutsche Bibliothek 

Die Deutsche Bibliothek lists this publication in the Deutsche Nationalbibliografie; 

detailed bibliographic data is available in the Internet at <http://dnb.ddb.de>. 



CR Subject Classification (1998): C.3, C.l.m, F.3, D.2, F.1.2, J.2, 1.6 
ISSN 0302-9743 

ISBN 3-540-21259-0 Springer- Verlag Berlin Heidelberg New York 



This work is subject to copyright. All rights are reserved, whether the whole or part of the material is 
concerned, specifically the rights of translation, reprinting, re-use of illustrations, recitation, broadcasting, 
reproduction on microfilms or in any other way, and storage in data banks. Duplication of this publication 
or parts thereof is permitted only under the provisions of the German Copyright Law of September 9, 1965, 
in its current version, and permission for use must always be obtained from Springer- Verlag. Violations are 
liable for prosecution under the German Copyright Law. 

Springer- Verlag is a part of Springer Science+Business Media 

springeronline.com 

© Springer-Verlag Berlin Heidelberg 2004 
Printed in Germany 

Typesetting: Camera-ready by author, data conversion by PTP-Berlin, Protago-TeX-Production GmbH 
Printed on acid-free paper SPIN: 10993194 06/3142 5 4 3 2 1 0 




Preface 



This volume contains the proceedings of the 7tlr Workshop on Hybrid Systems: 
Computation and Control (HSCC 2004) held in Philadelphia, USA, from March 
25 to 27, 2004. The annual workshop on hybrid systems attracts researchers 
from academia and industry interested in modeling, analysis, and implementa- 
tion of dynamic and reactive systems involving both discrete and continuous 
behaviors. The previous workshops in the HSCC series were held in Berkeley, 
USA (1998), Nijmegen, The Netherlands (1999), Pittsburgh, USA (2000), Rome, 
Italy (2001), Palo Alto, USA (2002), and Prague, Czech Republic (2003). This 
year’s HSCC was organized in cooperation with ACM SIGBED (Special Interest 
Group on Embedded Systems) and was technically co-sponsored by the IEEE 
Control Systems Society. 

The program consisted of 4 invited talks and 43 regular papers selected from 
117 regular submissions. The program covered topics such as tools for analysis 
and verification, control and optimization, modeling, and engineering applicati- 
ons, as in past years, and emerging directions in programming language support 
and implementation. The program also contained one special session focusing on 
the interplay between biomolecular networks, systems biology, formal methods, 
and the control of hybrid systems. 

We would like to thank the program committee members and reviewers for 
an excellent job of evaluating the submissions and participating in the on-line 
program committee discussions. Special thanks go to Edmund M. Clarke (Car- 
negie Mellon University, USA), John Doyle (California Institute of Technology, 
USA), Patrick Lincoln (SRI, USA), and Harvey Rubin (University of Pennsyl- 
vania, USA) for their participation as invited speakers. We are also grateful to 
the Steering Committee for helpful guidance and support. Many other people 
worked hard to make HSCC 2004 a success. We would like to thank T. John Koo 
for publicity, Janean Williams for local arrangements, and Valya Sokolsky for 
putting together the proceedings. We would like to express our gratitude to the 
US National Science Foundation and the University of Pennsylvania for finan- 
cial support. Their support helped us to reduce the registration fee for graduate 
students. 
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Abstract. We introduce the class of lazy rectangular hybrid automata. 
The key feature of this class is that both the observation of the continuous 
state and the rate changes associated with mode switchings take place 
with bounded delays. We show that the discrete time dynamics of this 
class of automata can be effectively analyzed without requiring resetting 
of the continuous variables during mode changes. 



1 Introduction 

We introduce here a class of linear rectangular hybrid automata called lazy hybrid 
automata and study its discrete time behavior. A central feature of this class is 
that the sensors report the current values of the variables and the actuators 
effect changes in the rates of evolution of the variables with bounded delays. 
More specifically, the state observed at X \ is a state that held at some time in 
a bounded interval contained in (Tfc_i, T*,). Further, if an instantaneous mode 
change has taken place at T then any necessary change in the rate of a variable 
will not kick in immediately. Rather, it will do so at some time in a bounded 
interval contained in (Tfc, Tk+i). A final restriction is that each variable’s allowed 
range of values is bounded. For convenience, we study the case where there is 
a single rate vector associated with each control state instead of a bounded 
rectangular region of vectors as is customary for rectangular hybrid automata [2] . 

Since both sensors and actuators have delays associated with them, a single 
symbolic trajectory of the automaton can correspond to uncountably many con- 
crete trajectories; even in a discrete time setting with the initial region being a 
singleton. Hence computing the discrete time behavior of a lazy hybrid automa- 
ton is non-trivial. Our main result is that this can be carried out effectively. It 
then follows that the discrete time behavior of a network of lazy hybrid automata 
that communicate by synchronizing on common actions can also be effectively 
computed. 

As is well known, the continuous variables available to an hybrid automaton 
and the fact that their rates of evolution can change instantaneously during a 
mode switch endows them with a great deal of expressive power. As a result, in a 
variety of settings, the control state reachability problem becomes undecidable, 
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as reported for instance, in [3]. A sharp characterization of the boundary between 
decidable and undecidable features of hybrid automata is provided in [8] as well 
as [2]. These results, as also the positive results reported elsewhere - for example, 
[5,13,11,10] - make it clear that except under very restrictive settings, one can 
not expect to get decidability if the continuous variables don’t get reset during 
mode changes; particularly in case their rates change as a result of the mode 
change. Viewed as a model of digital control systems that interact with physical 
plants through sensors and actuators, the resetting requirement severely restricts 
the modeling power of the automaton. Our results show that by introducing 
bounded delays into the functioning of the sensors and actuators, we can allow 
the variables to retain their values during mode changes. Admittedly, our positive 
results are obtained in the restricted setting of rectangular hybrid automata but 
the wealth of research concerning this setting (for instance, [6, 8, 5, 7]) suggests 
that this is a natural and well motivated starting point. 

We study the discrete time semantics of lazy hybrid automata. From a tech- 
nical standpoint, our work is a generalization of [7] where the discrete time 
behavior of rectangular hybrid automata is studied with the requirement that 
all instantaneous transitions should take place only at integer-valued instances 
of time. In our terms, [7] further assumes that the sensors and actuators function 
with zero delays which simplifies their analysis problem. In our setting, things 
are more complicated due to the non-zero delays associated with the sensing of 
values and actuating rate changes. Further, we feel that the approach proposed 
here is of some independent interest from a modeling point of view. It may also 
lead to the tractable analysis of larger classes of hybrid automata. Finally, our 
focus on discrete time semantics is relevant -as also argued in [7]- in that, as a 
model of digital controllers for continuous plants, the discrete time semantics of 
hybrid automata is more natural and useful than the continuous time semantics. 

Our work is, at least conceptually, in line with previous attempts to reduce 
the expressive power of timed and rectangular automata by taking away their 
ability to define trajectories with infinite precision [4,9,12]. Typically one de- 
mands the set of admitted trajectories to be “fuzzy”; if a trajectory is admitted 
by the automaton then it should also admit trajectories that are sufficiently close 
to the trajectory where “closeness” is captured using a natural topology over the 
trajectories. Surprisingly enough, this idea does not lead to more tractability as 
detailed in [9] and [12]. The key difference between our work and these previ- 
ous works is that in lazy hybrid automata, the fuzziness is introduced into the 
dynamics ; the observed continuous state based on which a mode change takes 
place at an instant is different from the actual continuous state that holds at 
that instant. Similarly, the actual rate at which a variable may be evolving at 
an instant may be different from the rate demanded by the current mode of the 
automaton. 

In the next section, we formulate the model of lazy hybrid automata. In 
section 3 we prove our main result, namely, the language of state sequences and 
action sequences generated by a lazy hybrid automaton are regular. Moreover, 
finite state automata representing these languages can be effectively computed. 
In section 4 we discuss the restrictions placed on lazy automata and point out 
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that many of them can be easily relaxed. We also sketch how our main result can 
be easily extended to networks of automata which communicate by performing 
common actions together. In the concluding section we summarize and point to 
some possible future work. 

2 Lazy Hybrid Automata 

Fix a positive integer n and one function symbol Xi for each i in {1,2,..., n}. 
We will view each as a function Xi : IR>o K >• 1R with 1R being the set of reals 
and IR>o, the set of non-negative reals. We let Q denote the set of rationals and 
I denote the set of closed intervals of the form [l, r] with l, r £ Q and l < r. We 
view [l, r] as the subset of 1R given by {z \ l < z < r}. 

A lazy hybrid automaton is a structure 
A = (Q, Act,q in ,V in , D,{p q } qeQ ,B, — >) where: 

— <5 is a finite set of control states. 

— Act is a finite set of actions. 

— tin € Q is the initial control state. 

— Vin G Q" is the initial valuation. 

— D = {g, 8 g , h, d/j} C Q is the set of delay parameters such that 

0<g<g + Sg<h<h + Sh<l- 

— p q € Q n is a rate vector which specifies the rate p q (i) at which each Xi 
evolves when the system is in the control state q. 

— B = B max ] £ I is the allowed range. 

— — Q x Act x I" x Q is a transition relation such that q ^ q' for every 
(q, a, /, q') in — >. Furthermore, if (q, a, /, q'), (q, a, I’, q') £ — > then I = I'. 

We shall study the discrete time behavior of our automata. At each time 
instant 71- , the automaton receives a measurement regarding the current values 
of the xfs. However, the value of Xi that is observed at Tk is the value that held 
at some t £ [Tk - 1 + h, Tk - 1 + h + Sh\- If at Tk, the automaton is in control state 
q and observed n-tuple of values (fi, ^ 2 , • • • , v n ) is in / with (q, a , /, q') being a 
transition, then the automaton may perform this transition instantaneously by 
executing the action a and move to the control state q' . Thus as usual, the xfs 
will cease to evolve at the rates p q and instead start evolving at the rates p q >. 
However, this change in the rate of evolution will not kick in at Tk but at some 
time t £ [T k + g,T k + g + <$ g ]. In this sense, both the sensing of the values of 
the xfs and the rate changes associated with mode switching take place in a 
lazy fashion but with bounded delays.. We expect g to be close to 0, h to be 
close to 1 and both 8 g and Sh to be small compared to 1 so that in the idealized 
setting, the change in rates due to mode switching would kick in immediately 
(g = 0 = S g ) and the value observed at Tk is the value that holds at exactly Tk 
[h= 1 and Sh = 0). Indeed, this is the setting considered in [7]. 

B specifies the range of values within which the automaton’s dynamics are 
valid. The automaton gets stuck if any of the xfs ever assume a value outside the 
allowed range [B m i n , B max \. A number of the restrictions that we have imposed 
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are mainly for ease of presentation. Later, we will discuss how these restrictions 
can be relaxed. Our main result is that the control state and action sequence 
languages generated by a lazy hybrid automaton are both regular. Furthermore, 
these language can be computed effectively. 

2.1 The Transition System Semantics 

Through the rest of this section we fix a lazy hybrid automaton A as defined 
above and assume its associated notations and terminology. The behavior of A 
will be defined in terms of an associated transition system. 

A valuation is just a member of 1R™. We let i range over {1, 2, . . . , n}. The 
valuation V will be viewed as prescribing the value V(i) to each variable Xi . 
A configuration is a triple (g, V, q') where q, q' are control states and V is a 
valuation, q is the control state holding at the current time instant and q' is the 
control state that held at the previous time instant. V captures the actual values 
of the variables at the current instance. The configuration (g, V, q') is feasible iff 
V (l) ^ [Bmim B max ] for every i. The initial configuration is, by convention, 
the triple (qimVimQin)- We assume without loss of generality that the initial 
configuration is feasible. We let C a denote the set of configurations. Since A 
will be clear from the context, we will often write C instead of C 4 . 

As in the case of timed automata [1], a convenient way to understand the 
dynamics is to break up each move of the automaton into a time-passage move 
followed by an instantaneous transition. At To, the automaton will be in its 
initial configuration. Suppose the automaton is in the configuration ( Qk ■ Vk , q' k ) 
at Tfc. Then one unit of time will pass 1 and at time instant T k+ 1 , the automaton 
will make an instantaneous move by performing an action a or the silent action 
r and move to a configuration (q k +i, Vfc+i, q' k . x ). The silent action will be used 
to record that no mode change has taken place during this move. Again, as 
often done in the case of timed automata, we will collapse the two sub-steps 
of a move (unit-time-passage followed by an instantaneous transition) into one 
“time-abstract” transition labeled by a member of Act or by r. 

With this intuition in mind, we now define the transition relation 
=>C C x Act U {r} x C as follows. 

— Let (q,V,q r ) and (gl,Vl,gl') be configurations and a £ Act. Then 
(q,V,q r ) ==> {ql,Vl,ql') iff qY = q and there exists in A a transition 

of the form q -^4 ql and there exist tl £ [g,g + S g ] n and t2 £ [h, h + Sh] n 
such that the following conditions are satisfied for each i. 

(1) VI (i) = V(i) +j» q '{i) ■ tl(i) +J> q (i)-( 1 ~ 

(2) V(i) + Pq'{i ) • tl(i) + p q (i) ■ (t2(i) — fl(i)) £ I(i ) for each i. 

— Let {q,V,q r ) and (gl, VI, ql') be configurations. Then (q, V, q’) =A> 
(gl, VI, ql') iff gl = ql' = q and there exists tl £ [g, g + 8 g ] n such that 

VI (i) = V{i) + p q '(i ) • tl(i) + p q {i) ■ (1 — £l(i)) for each i. 

1 We assume that the unit of time has been fixed at some suitable level of granularity 
and that the rate vectors {pq\ have been scaled accordingly. 
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Basically there are four possible transition types depending on whether q = q' 
and a £ Act. Suppose (q,V,q') ==> (ql,Vl,ql') with a £ Act. Assume that 

q —4 ql in A and tl £ \g, g + 6 g ] n and t,2 £ [h, h + Sh] n are as specified above. 
We first note that ql ^ q by the definition of the transition relation of A. The 
requirement ql ' = q captures follows from our convention that ql' is the control 
state that held in the previous instant and we know this was q. 

First consider the case q A q' and let us suppose that the configuration 
(q,V,q') holds at T k . We take q ^ q' to mean that a change of mode from 
q' to q has just taken place (instantaneously) at T k based on the observations 
that were made available at T k . However, at T k , the automaton will continue 
to evolve at the rate dictated by p' q . Indeed, each Xi will, starting from T k , 
evolve at rate p' q (i) until some Tk + t\ with t\ £ \g, g + S fl }. It will then start 
to evolve at rate p q (i) until T k +\. Consequently, at T k + 1, the value of will be 
Vl(i) = V(i) + p' q (i) ■ ti + p q (i) ■ (1 — H). On the other hand, ql ^ q implies 
that another instantaneous mode change has taken place at T k + 1 based on the 
measurements made in the interval [T k + h, T k +h+Sh\- Suppose Xi was measured 

at Tk+t2 with f 2 £ [Tk + h, T k + h+6h\. Then in order for the transition q -^4 ql 
of A to be enabled at T k + 1, it must be the case that the observed value of a 
at Tk + t'2 falls in I(i). But then this value is V(i) + p' q (i) ■ t\ + p q (i) ■ (0 — t\). 
This explains the demands placed on the transition (q,V,q') (ql, VI, ql'). 
It is worth noting that in case q = q' (i.e. no mode change has taken place at 
Tk) then VI (i) = V(i) + p q (i) ■ h + p q (i ) • (1 — H) = V (i) + p q as it should be. 
Furthermore, V(i) + p q (i) ■ t\ + p q (i) • (t 2 — H) = V(i) + p q (i) ■ t 2 and this must 
fall in I(i) as to be expected. 

Similar (and simpler) considerations motivate the demands placed on transi- 
tions of the form (q, V, q') =^=> (ql,Vl, ql'). Here again, it is worth noting that, 
in case q = q' , VI (i) is determined uniquely, namely, VI (i) = V(i) + p q (i). 

We now define the transition system 
TSa — (AC(4, (qin, Vi n , qin) , Act U {t}, — s ^4) via. 

— RCA, the set of reachable configurations of A is the least subset of C that 

contains the initial configuration (f/,;„, V in , q ln ) and satisfies: 

Suppose ( q , V, q') is in RC_ 4 and is a feasible configuration. Suppose further, 

(q, V, q') (ql, VI, q) for some a £ Act U {r}. Then (ql,Vl, q) £ 

— =>A is => restricted to RCa x Act U {r} x RC 

We will often write RC instead of RC, 4 and write => instead of =>a- We 
note that a reachable configuration can be the source of a transition in TSy\ 
only if it is feasible. Thus infeasible reachable configurations will be deadlocked 
in TSa- 

A run of T Sa is a finite sequence of the form 
o' = (qo, V 0 , q' 0 ) a 0 (qi, Hi, q[) ai (q2,V 2 , q' 2 ) ■ ■ ■ (%, V k , q' k ) where (q 0 , V 0 , q' 0 ) is the 
initial configuration and (qm,Vm,q' m ) (q m +i, V m+1 , q' m+1 ) for 0 < m < k. 
The st-sequence (state sequence) induced by the run 0 above is denoted as st(a) 
and it is the the sequence qoqi ■ ■ ■ q n . On the other hand, the act-sequence induced 
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by a is denoted as act (a) and it is the sequence aooii . . . a n . We now define the 
languages C st (A ) and C act (A) as : 

— C st (A) = (st(cr) | a is a run of A}. 

— C ac t{A) = {act(cr) | a is a run of A}. 

Our main result is that jO st (A) is a regular subset of Q* while C ac t{A ) is a reg- 
ular subset of ( Act U {r})*. Moreover, we can effectively construct finite state 
automata representing these languages. As a consequence, a variety of verifica- 
tion problems for lazy hybrid automata can be effectively solved. 

3 Proof of the Main Result 

We shall first establish the main result for the one dimensional case. As is often 
the case with rectangular hybrid automata [5], it will then be easy to lift the 
proof to the n-dimensional case with the help of a (Cartesian) product operation. 



3.1 The One Dimensional Case 

Let A = (Q, Act,qi n ,Vi n , D, {p q } q& Q,B , — >) be a lazy automaton. We assume 
for A , the terminology and notations defined in the previous section. Until fur- 
ther notice , we set n = 1 and we will write x instead of Xi and p q instead of p q {i) 
for q € Q. The key idea is quantize the unit time interval and correspondingly 
the phase interval [ B min , B max ], We first define A to be the largest positive 
rational number that integrally divides every number in the finite set of rational 
numbers {g,8 g ,h,8h,l}- We next define T to be the largest positive rational 
number that integrally divides each number in the finite set of rational numbers 
{ Pq * A | q £ U Bmax } U {/, r | ( q , a , [l, r\,q') is a transition in A}. 

Let 7Z denote the set of integers. We now define the map 
||-|| : IR — > 7Z x ({0, 1} U {_L}) as follows. 

- If V £ (- 00 , B min ), then ||u|| = ( k min - 1, _L) where k min ■ T 
v S (B max , oo) then ||u|| = (k ma xA L) where k max ■ T = B max . 

— Suppose v € [B m i n , B max \ and k € 7Z and v € [0, T) such that v 
Then ||v|| = (k, 0) if v = 0 and ||u|| = ( k , 1) if v ^ 0. 

The map ||-|| can be extended in a natural way to configurations. Denot- 
ing this extension also as ||-||, we define ||(g, v, < 7 ') || to be (q, ||u|| , q'). Let 
Ba = { ||c|| | c £ Ca }• Clearly T>a is a finite set and we will often write V 
instead of T>a- Our goal is to show that the equivalence relation over the reach- 
able configurations RC of A induced by the map ||-|| in turn induces a right 
congruence of finite index over Q*. The proof will make use of the following 
simple observation. In stating the observation and elsewhere, we will use the 
following notations. For q,q' £ Q we let N q and N qq > be the positive integers 
such that | p q ■ A\ = N q ■ r and | (p q — p' q ) ■ A\ = N qq f ■ T. Clearly, N q and N qq i 
exist because of the choice of A and r. 



Bmin- 1 ^ 
= k- T +V. 
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Lemma 1. Let q,q' £ Q. Define the functions f q and f qq ’ as: 

(1) fq ■ [0, A/N q ] [0, r] and is given by f q (9) = \p q ■ 9\. 

(2) .fqq' ■ [0, A/N qq /] [0, r] and is given by f qq '{9) = | (p q - p q >) ■ 9 1. 

Then both f q and f qq > are well-defined, continuous and onto. 

Proof. Follows easily from the definitions and the basic property of monotonic 
real valued functions over bounded domains. □ 

We are now ready to tackle the main part of the proof. 

Theorem 1. Let cl and c2 be two reachable configurations such that ||cl|| = 
1 1 c2 1 1 . Suppose a £ Act U {r} and cl' is a reachable configuration such that 
cl ==>a cl'. Then there exists a reachable configuration c2' such that c2 ==>.4 c2' 
and Hcl'll = ||c2'|| . 

Proof. Clearly cl is feasible and since ||cl|| = ||c2||, it follows that c2 is also 
feasible. 

Assume that cl = {ql,Vl,ql') and c2 = (q2,V2,q2') and that ||V1|| = 
(ATI, zl) and ||V2|| = ( K2,z2 ). Since ||cl|| = ||c2||, we can set q = ql = q2, 
q' = ql' = q2! and ( K,z ) = (Kl,zl) = ( K2,z2 ). If 2 = 0 then VI = V2 and 
hence cl = c2 and the result follows. 

So assume that z = 1 and VI ^ V2. Hence VI, V2 £ ( K.r,(K + 1 ).r) 
and hence || (<7, HI, g')|| = ||(g, V2, q')\\ = (q,(K,l),q'). Furthermore, there exist 
vl, v2 £ (0, r ) such that vl ^ v2 and VI = K ■ P + vl and V2 = K ■ T + v2. 
In what follows, for the sake of convenience, we will assume that 0 < p q ’ < p q 
and that v2 < vl. From the structure of the proof it will be obvious that this 
involves no loss of generality. 

Let cl' = (s,VT, q). Then we have (q, Vl, q') ==> (s, Vl' , q). We are required 
to show that there exists V2' such that {q,V2,q') (s,V2’,q) with || Vl 7 1| = 

|| V2 ' || . We shall do this by considering four cases. 

Case 1: q = q' and a = r. 

Since q = q ' , no mode change has taken place in the previous time interval. 
Hence the automaton will evolve at rate p q during the current unit interval. On 
the other hand, a = r implies that s = q and hence no mode change takes place 
at the end of this unit interval either. Consequently, we must have Vl' = Vl+p q . 
We now set V2' = V2 + p q . Then it follows that {q,V2,q') (q,V2',q). We 

need to argue that ||V1 , || = ||V2'||. 

In what follows, we define for £ £ {g, S g , h, 5h, 1}, to be the positive integer 
satisfying £ = ■ A. These positive integers must exist by the choice of A. 

Now p q = p q ■ 1 = p q ■ Ni ■ A = N q ■ Ni ■ r. (Recall that p q ■ A = N q ■ T). But 
then Vl, V2 £ (K • T, (K + 1) • T) and hence Vl', V2' £ ((K + N q ■ N x ) ■ T, (K + 
1 + N q ■ Ni) ■ r). This at once leads to ||V1'|| = ||V2 , ||. 

Case 2: q = q' and a £ Act. 

Since q = q' we again have that no mode change has taken place in the 
previous interval and hence the automaton will evolve at rate p q in the current 
interval. Hence, as in the previous case, we must have VV = Vl + p q . Again, we 
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set V2' = V2 + p q . Consequently as shown in the previous case, ||V1'|| = ||V2'||. 
So if we show that (q, V2, q') =^> (s, V2', q ), then we are done. 

We are given that (g, Vl, q') (s, VI', q). Hence there exists a transition of 

the form (g, a, I, s) in A and there exists 1 1 £ [h, h+Sh] such that Vl+p q -tl £ I. 
We just need to show that there exists t2 £ [h, h + 8h\ such that V2 + p q -t2 £ I. 

In order to fix t2, recall that h = Nh ■ A and 8h = Ng h ■ A. We first note that 
tl £ [N h -A, (Nh+Ng h )-A\. Noticing that p q -A = N q -T and hence p q -{A/N q ) = 
r we set A q = A/N q , and observe that tl £ [Nh ■ N q - A q , ( Nh + N$ h ) ■ N q ■ A q \. 
Let N be the least integer in the interval [Nh ■ N q , (N h + Ng h ) ■ N q ] such that 
tl £ [N ■ A q , (N + 1) • A q }. Let 61=tl-N- A q . Clearly 01 £ [0, A q [. 

Suppose 01 = 0. Then p q ■ tl = p q ■ N ■ A q = N ■ T and hence VI = 
Vl + p q -tl £ ((K + N)-r, (K + l + N)-T). Set t2 = tl. Then V2 = V2 + p q -tl £ 
{{K + N) ■ r, ( K + 1 + N) ■ r) too. Now assume that / = [l, r\. Then there exist 
integers Ni and and N r such that l = Ni- T and r = N r ■ T with Ni < N r . Since 
VI £ [l, r], we must have N < ( K + N) < ( K + N + 1) < N r . But this implies 
that V2 = V2 + p q ■ tl £ [l, r] too. Hence (q,V 2, q') =^> (s, V2', q). 

The case 01 = A q can be dealt with in a similar manner by again setting 
t2 = tl. 

So now assume that 01 £ (0, A/N q ). Then clearly VI = VI + p q ■ tl £ 
[vl + (K + N)-r,vl + (K + N+l)-r]. (Recall that vl = VI — K ■ T and 
v2 = V2 — K ■ r.) There are three possibilities to consider. 

Firstly, suppose VI £ [vl + ( K + N) ■ T, (K + N + 1) • T). Then we set 
t2 = N-A q . Clearly V2 = V2 +p q -N-A q £ ((K+N)-r, (K+N+l)-r). But then 
Vl £ [vl + (K+N)-T, (K+N+l)-r) implies Vl £ ((K+N)-r, (K+N+l)-r). 
Consequently Vl £ [l, r] implies N[ < ( K + N) < [K + N + 1) < N r as before 
and this in turn implies V2 £ [l, r]. This leads to (q,V 2, q') (s, V2', q). 

Secondly, suppose vl = (K + N + 1) • V. Then, ( K + N + 1) • V £ (v2 + ( K + 
N) ■ r, v2 + ( K + N + 1) • r). From Lemma 1, it follows that there exists 02 in 
[0, Aq] such that v2 + {K + N)-r + p q - 0 2 = (K + N+l)-T. Set t2 = N ■ A q + 0 2 . 
Clearly, V2 = V2 + p q ■ t2 = Vl = (K + N + 1) • F. Again, Vl £ [l, r] implies 
V2 £ [l, ?’] as required. 

Thirdly, suppose Vl £ ((K + N + 1) • V, vl + ( K + N + 1) ■ T], Then we 
set t2 = {N + 1) • A q . Clearly V2 = V2 + p q ■ {N + 1) • A q £ {(K + N + 1) • 
r, (K + N + 2) ■ r). But then Vl £ (( K + N + 1) ■ T, vl + {K + N + 1) ■ T] 
implies vl £ (( K + N + 1) • V, (K + N + 2) • V). Thus again, Vl £ [l, r] implies 
V2 £ [l, r]. 

Case 3 ■■ q^ q' and a = r. 

Since q ^ q ' , an instantaneous transition has taken place at the end of the 
time-passage move leading to (q,Vl,q r ). Hence the automaton will continue to 
evolve at rate p q ' until some tl £ [g, g+8 g ] and then will evolve at the rate p q for 
the rest of the period 1— tl. Moreover tl is such that Vl' = Vl+p' q -tl+p q -(l— tl). 
We need to find t2 £ [g, g + <5 S ] such that V2' = V2 + p q > ■ t2 + p q ■ (1 — t2) and 
|| Vl'H = || V2'|| . In order to fix t2, let g = N g ■ A and S g = Ns g ■ A. 
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Noticing that (p q — p q ') ■ A = N qq / • A and hence (p q — p q >) ■ ( A/N qq /) = A we 
set A qq , = A/N qq i, and observe that £1 £ [AT a • A^/ • A qq > , (A^ + A^J-A/^/- A ? g']- 
Let AT be the least integer in the interval [N g ■ N qq >, (N g + Ng ) ■ N qq >] such that 
£1 £ [N-A qq ., (N + 1) • A qq >]. Let 91 = tl — N ■ A qq >. Clearly 91 £ [0, A,,-]. 

We now have VI' = VI + p q > ■ N ■ A qq > + p q > ■ 91 + p q ■ (A qq > - 91) + p q ■ (Ni ■ 
N qq i — N — 1) • A qq i . (Recall that N\ ■ A = 1.) Expanding this expression and 
simplifying using the definitions of N q , N q i, N qq t and A qq i, we get: 

VI' = Vl+(N v N q -N)-r-(p q -p q ,)-91. We recall that vl = Vl-K-Tand v2 = 
V2-K- r. Since 91 ranges over [0, A qq >], we have that (p q — p q i) ■ 91 ranges over 
[0, r). Hence VI' £ [vl + (K+N!-N q -N)r, vl + (K+N!-N q -N+i)-r]. Again 
there are three situations to consider. For convenience, let K' = N\ ■ N q — N . 

Suppose V 1' £ [vl + (AT + K') • T, (K + K' + 1) • A) . Then we set t2 = N ■ A qq > . 
Then it is easy to see that t2 £ [g,g + S q \. Now let V2! = V2 + p q > -t2 + p q - (l — t.2). 
Then by our choice of t2, we have, V21 = V2+p q fN-A qq '+p q -(Ni -N qq '— N)-A qq i. 
Simplifying this expression, we get V2' = V2 + K' ■ A. Since V2 = v2 + K ■ A, 
we then get V21 £ ((AT + AT') • A, ( K + K' + 1) • A). As a result, ||V1'|| = ||V r 2'||. 
By the choice of t2, it is also clear that (q,V2,q r ) (s,V2',q). 

The case El' £ ((AT + K' + 1) • A, vl + (k + AT + 1) • A] is handled in a similar 
manner by setting £2 = ( AT + 1) • A qq r . 

So suppose that VV = ( K + K' + 1) • A. Then by Lemma 1 we can find 
92 £ (0, A qq i ) such that with £2 = N ■ A qq i + 92 , and V2' = V2 + p q > ■ £2 + 
p q • (1 — £2), we can obtain V2' = ( K + K' + 1) • A. This follows from the 
fact that as 92 ranges over [0, A qq i], we will have V2! ranging continuously over 
[v2+(A' + AT , )-A, v2 + (A' + A' , + l)-A] and surely (K + K' + 1)- A lies within this 
range. Clearly by the choice of £2 and V2', we have (q,V2,q r ) =?==> (s,V2',q). It 
also follows at once that ||C1 , || = ||C2 , ||. 

Case 4 : q^ q' and a £ Act. 

This is the most general case where the rate will change during the current 
period and the time-pass move will be followed by an instantaneous execution 
of a transition of A. 

Since ( q , VI, q') =A> (s, VV, q), there exist £1 £ [ g , g+S g ] and £1' £ [h, h+Sh) 

and a transition q 9 s in A such that VV = VI + p q > • tl + (1 — £1) • p q 

and VI + p q i ■ tl + p q ■ (tl' — tl) £ I. We need to find £2 £ \g, g + 5 g \ and 

£2' £ [h, h + Sh] such that V2 + p q > ■ £2 + p q ■ (£2' — £2) £ / and ||d'|| = ||C2'|| 

where V2' = V2 + p q > • £2 + (1 — £2) • p q . 

As before, we set A gq i = A/N qq i and let N be the least integer in the 
interval [N g -N qq i, (N g + Ng g ) ■ N qq i] such that £1 £ [N-A qq i, (N + 1) ■ A qq i], Let 
91 = tl — N ■ A qq i . Clearly 91 £ [0, A qq i], Using the argument developed to settle 
the previous case, we can conclude that VI' = Vl+(Ni-N q —N)-r—(p q —p q i)-61. 
As before, we set K ' = N\ ■ N q — N. We need to examine two cases. (It is worth 
recalling here that we are operating under the assumptions 0 < p' q < p q and 
v2 < vl). 

Suppose VI' £ [vl + (K+K')-r, v2 + (K+K' + 1)- A]. Consider £2 = (AT+1)- 
A qq i +92 for some 92 £ [0, A qq ,]. Define V2' = v2+(K'+K+l)-r~92-(p q -p q i). 
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As 92 ranges over [0, A qq >], V2' will range over [v2 + ( K ' + K) ■ T, v2 + (I\ + 
K' + 1) • A]. Hence, by Lemma 1, we can fix a 92 such that V2' = HI'. 

Suppose on the other hand, VV £ (v2 + (K + K' + 1) ■ T, vl + (K+K' + 1)-T]. 
Then we set 92 = 0 so that t2 = ( N+1)-A qq > and hence V2' = v2+(K+K' +1)-T . 
Clearly both VV and V2! lie in (( K + K' + 1) • T, (K + K' + 2) ■ r). Hence 

Ill'll = ll^2'||. 

We note that in either case, our choice of 92 guarantees that VV = V2! or 
V2! < VV with VV - V2' <vl- v2. 

Turning to the choice of t2', we define as before, A q = A/N q . Let J be 
the least integer in the interval [Nh ■ N q , ( Nh + Ng h — 1) • 7V g ] such that t V £ 
[J-A q , (J+ 1) • A q \. Let 9V = tl' - ( J • A q ). Clearly 9V £ [0, A q \. 

Let HI" = Vl+p ql .tl+p q -(tl'-tl). Then HI" = Vl+p q '-N-A qq ,+p ql -91 + 
p q ■ (. A qq > - 91) + p q ■ ( J ■ N q ■ A q - (N + 1) • A qq i) + p q ■ 6V . Again expanding and 
simplifying this expression, we get HI" = HI + (N q ■ J — N) ■ P — (p q — p q f) ■ 91 + 
p q -9V. Let L = N q - J — N. Then HI" = vl+(K + L) ■ T - (p q - p q ,) -91 + p q -9V. 

Now HI"' = vl+(K+L)-r~ (p q — p q ’)-91 must lie in [vl+(K+L)-r, vl+(K+ 
L+l)--T]. Suppose HI'" lies in [vl + (k+L)-r, v2+(K+L+l)-r], Then our choice 
of 92 ensures that u2 + (A'+L)-P— (p q — p q ’)-92 = vl + (K + L) ■ T — (p q — p q r) ■ 91. 
We now set 92' = 9V and t2' = J ■ A q + 92! . Clearly t2' £ [h, h + 5f,] and 
H2" = H2 + p q i ■ t2 + p q ■ (t2' — t2) = HI" £ I and hence we have, as required, 
(q,V2,q') (s,V2’,q) with ||H1'|| = ||H2'||- 

Finally, assume that HI"' = vl + {K + L) ■ T — (p q — p q >) ■ 91 lies in (v2 + 
(K + L + 1) • r, vl + (K + L + 1) • r] Then our choice of 92 ensures that 
H2"' = v2+ (K + L)T - ( p q - Pq ,) • 92 = v2 + (K + L + 1) • T and thus 
HI"' — H2'" < vl — v2. Now depending on 9V, the value of HI" must lie in 
(v2 + (K + L + 1) • T, vl + (K + L + 2) ■ T]. If HI" lies in (v2 + (k + L + 1) • 
r, v2 + (k + L + 2) • T] then we can, by lemma 1, pick 92! £ [0, A q \ so that 
H2" = HI" where H2" = H2 + p q , ■ t2 + p q ■ (t2' - t2) with t2' = J ■ A q + 92'. If 
on the other hand, HI" lies in (v2 + ( I\ + L + 2) ■ T, vl + (k + L + 2) ■ T] we can 
set 92! = A q and t2' = J ■ A + 92' so that H2" = v2 + ( K + L + 2) ■ T. In either 
case, we have t2' £ [h, h + Sh] and H2" £ / so that (q,V2,q') (s, V2',q) 

with || HI' || = ||H2'||. □ 

We now define the finite state automaton Zj± = (T>, (qi n , (ko, 0), qi n ), Act U 
{r}, V) where ko-T = H„ and the transition relation T>x (ActU{r}) xT> is 
given by: (q, ( k , d),q') (ql, (. kl , dl),ql') iff there exist configurations (q, H, q') 

and (gl,Hl,gl') such that (q,V,q') (gl,Hl,gl') and ||H|| = ( k,d ) and 
|| HI || = (kl, dl). In what follows, we will often write as just Z. Note that, we 
are setting all the states of Z to be its final states. 

We define C st (Z) to be the subset of Q* as follows. A run of Z is a se- 
quence of the form (q 0 , (l 0 ,d 0 ),q' 0 ) a 0 (qi, (h, di), q[) aq . . . (q m ,(l m ,d m ),q' m ) 
where (q 0 ,(l 0 ,d 0 ),q' 0 ) = (q in , (k 0 , 0), q in ) and 

(qj, (lj, dj), q!) (qj+i, (lj+i, dj+i), q' j+1 ) for 0 < j < m. Next we define 
qodi ■ ■ -q-m G Cst(Z) iff there exists a run of Z of the form 

(q 0 ,(l 0 ,d 0 ),q' 0 ) a 0 (qi, (h, di), q[) a x ... (q m ,(l m ,d m ),q' m ). Clearly C st (Z) is a 




Lazy Rectangular Hybrid Automata 



11 



regular subset of Q* and it does not involve any loss of generality to view Za 
itself as a representation of this regular language. 

Theorem 2. The automaton Z _ 4 can be computed effectively. Moreover 
C st (A) = C s t(Z^) and C ac t(A) = C(Z_ 4) where C(Za) is the regular subset 
of ( Act U {t})* accepted by Z a in the usual sense. (Note that all the states of 
Zjx are final states.) 

Proof. Clearly the finite set of states T> and the initial state ( q in , (k 0 , 0),q in ) can 
be computed easily. The transition relation -w is expressible in the first order 
theory of the real ordered field which is a decidable theory [14] . 2 . For instance, to 
determine if ( q , (k, 1 ), g') (gl, (kl, 1), g), with a £ Act, we first check if there is 
a transition tr in A of the form (g, a, I, ql). If there is no such transition then we 
conclude that ( q , (k, 1 ),q') (ql, (kl,dl),q) is not a transition in Za- If there 
is such a transition then for each such transition tr we construct the formula 
tptr, take the disjunction of all such formulas and check for its satisfiability. 

Suppose tr = (q, a, I, ql). Then ip tr will conjunctively assert the following: 

— There exists V such that k ■ T <V < (k + 1) • T. 

— There exists tl such that g < tl < g+S g and kl-T < V+p q >-tl+p q -(l — tl) < 

(ki + i)-r. 

— There exists t2 such that h < t2 < h+Sh and l < V + p q ' -tl + p q -(t2 — tl) < r 

(where I = [l, r]). 

To see that C st (A) = C st (Z) we hrst note that C st (A) C C st (Z ) follows from the 
definition of Z_ 4. To conclude inclusion in the other direction, we will argue that 
for each run (q 0 , (l 0 , d 0 ),q' 0 ) a 0 (qi, (h, di), q[) an . . . (q m , (l m , d m ),q' m ) of Z there 
exist V 0 , Vi . . . V m £ 1R such that (q 0 , V 0 ,q' 0 ) a 0 (qi,Vi,q{) ai . . . (q m , V m , q ’ m ) is 
a run of TSa- And furthermore, ||Vjj| = ( lj,dj ) for 0 < j < m. The required 
inclusion will then follow at once. For m = 1, it is clear from the definitions and 
so suppose that (q 0 , (l 0 ,do),q' 0 ) a 0 (q 1, (Zi,di),gj) . . . (q m , (l m ,d m ),q' m ) 
a m (qm+i,(lm+i,d m+ i),q' m+ i) is a run of Z. By the induction hypothesis, 
there exists a run (q 0 , V 0 , q' 0 ) a 0 (qi, V\, gj) ati . . . (q m , V m , q ' m ) of TSa with the 
property, ||Vj|| = ( lj,dj ) for 0 < j < m. 

Now (g m , (lm, d m ), q'm) ^ (g m+ i, (l m+ i,dm+i),q' m+ i) implies that there ex- 
ist TA and V^ +1 such that (g m ,V^,g^) ^ (q m+1 ,Vf n+1 ,q , m+1 ) and ||V^|| = 
(l m , d m ) and || V^ +i || = (l m +i, d m+1 ). But this implies that ||F4|| = ||F m ||. Hence 
by Theorem 1, there exists V^+i such that (qm,V m ,q' m ) ^ (q m+1 ,V m+1 , q' m+1 ) 
and moreover ||V^ l+1 || = ||F m+ i||. Thus C st (A) = C st (ZA). It now also follows 
easily that C ac t(A) = C(Zjf). □ 

In what follows, we will refer to Z as the zone version of A. 

2 This is an overkill as detailed later 
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3.2 The n-Dimensional Case 

We now consider an n-dimensional hybrid automaton A defined as in the pre- 
vious section with the associated terminology and notations. Our goal is to 
show that C st (A) is a regular subset of Q* while C ac t(A) is a regular subset of 
( Act U {r})*. 

To do so, we first define the family of one dimensional automata {-4*} = 
(Q, Act, q\ n , V? n , D, {p\} q&Q , B, — »*) where: 

— Wn(*) is the i - th component of Vi n . 

~ Pq= Pq{i) 

— q (-—>1 q' iff there exists q ^—5 q' in A with P = I(i). Again, I(i ) denotes 
the i-tlr component of I. 

Let Z l be the zone version of A‘ with 

Z l = (V 1 , (q in , (k l o ,0),q in ), Act U {r},~*j). We now define the finite state au- 
tomaton Z A = ( V , (qin, no, qin), Act U {r}, V) which will constitute the zone 

version of the n-dimensional automaton A as follows. 

— T>, the states of this automaton, will be of the form ( q,K,q ' ) with q,q' £ 
Q and k £ ((2Z x {0,1})". Let k = ((hi, d\), (& 2 , d^) ■ . ■ , (k n , d n )). Then 
(q, k, q') £ V iff there ( q , ( fc* , di),q') £ V 1 for each i in {1,2,..., ?r}. 

~ k 0 = ((fco,0),(A:o,0),...,(feJ,0) 

— V x ( Act U {t}) x V is given by: 

Let (q,K,q'),(ql,Kl,qV) £ T> with k = ((fci, di), (^ 2 , ^ 2 ) • • • , (k n , d n )) and 
nl = ((fell, dli), (fcl 2 ,dl 2 ) •••, (fel„,dl n )). Then (q,K,q') ((ql,Kl,ql') 
iff (q, (hi, di), q') (ql, (kli,dli),q) for each i £ { 1 , 2 , . . . ,n}. 

As before, we will often write Z instead of Z . 4 and refer to it as the zone 
version of A. We denote by C s t(Z) the state sequence language of Z and define 
it in the obvious way. We also define C(Z) to be the subset of (Act U {r})* 
accepted by the finite state automaton Z. 

Theorem 3. The automaton Z 4 can be computed effectively. Moreover 
C st (A) = C st (Z A ) and C ac t(A) = C(Z A ). 

Proof. Since, by Theorem 2, each of the finite state automata Z l can be com- 
puted effectively, so can Z be. The proof of the facts C st (A) = C st (Z A ) and 
Zact(A) = C(Z A ) is routine and we omit the details. □ 

As for the complexity of our decision procedure, we first estimate the size of the 
automaton and the time complexity of constructing the automaton for the one 
dimensional case. Let / be the total number of relevant intervals on 1R. In other 
words, I = (B max — B min )/r + 2. Then the number of states is 0(m 2 • I) where 
m = | Q | is the number of control states of the lazy automaton. For constructing 
the transitions, we need to check if there is a transition from (q, (k,l), q') to 
(<?1, (k 1, 1 ),q) labeled with the action a. It is clear that the most complex case 
is when a £ Act and q ^ q' and we need to check for the existence of at most 
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0(m 4 ■ 1 2 • \Act\) such possible transitions. To decide if such a transition exists 
for a given pair of states (q,(k, l),g') to (gl, (fcl, 1), q) and a given symbolic 
transition in the lazy automaton of the form (g,a, [l,r],ql) we need to check if 
there exists V and fl and t2 such that: 

- k ■ T < V < (k + 1) • r 

- g < tl < g + Sg 

- ki ■ r < v + /v • ti + p q • (l - ti) < (fci + 1 ) • r 

/ / A t2 A / 1 Sfi 

- I < V + pq’ ■ t\ + Pq ■ (t2 — tl) < r. 

The above are 10 linear inequalities in three variables V, tl, and 12. Linear 
programming allows us to check if they all can be satisfied. This can be done 
in time proportional to the length of the constraints as there are a constant 
number of variables and constraints. Therefore, the time to check each inequality 
is 0(log(j;)) since r requires the largest number of bits to represent of all the 
quantities in the inequalities. Thus the time complexity of building the automata 
in the one dimensional case is 0(m A ■ 1 2 • \Act\ • log(-p)). 

For the n dimensional case, the number of states is 0(m 2 ■ /”). So the num- 
ber of possible transitions is 0(m 4 • I 2n ■ |Act|). Now there are n groups of 10 
inequalities with each group involving 3 variables. Hence overall time complexity 
of constructing the automata is 0(m A ■ I 2n ■ \ Act\ ■ n ■ log( -p)). 

4 Some Extensions 

In order to simplify the initial presentation, we placed a number of restrictions 
on our automata. Here we first examine which of these can be relaxed so that, 
with minor overhead, our main results go through smoothly. We then formulate 
a composition operation for lazy hybrid automata in a standard way using which 
large automata can be presented in a succinct fashion. These networks of lazy 
hybrid automata can also be analyzed effectively. 

Let A = (Q, Act, qi n ,Vin, D, {pq} q( zQ, B, — ►) be a lazy hybrid automaton. 
We could permit a set of initial control states and a set of initial valuations 
for each initial control state, provided they can be specified using rectangular 
constraints. Our results will go through with minor modifications. It is also clear 
that our demand 0<g<g + 5 g <h<h + 5h<l is only for convenience. 
We could have different delay parameters for different variables and these delays 
could spill over more than one time unit. 

The restriction that there is at most one a-labeled transition between a pair 
of control states is mainly for convenience. If this condition is violated we could 
use renaming to enforce this property, construct the zone automaton and then 
restore the old names. 

State invariants can be introduced in the expected manner and we could al- 
low resets of the variables during a mode switch. Finally, we have avoided the 
customary use of differential inclusions to specify the rates mainly to avoid clut- 
ter. Our results will still go through, with some additional notational overhead, 
if we permit these extensions. 
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The boundedness restriction on the allowed range B = [B m i n , B max \ is crucial 
though, from a modeling point of view, it is not crippling. The fact that we have 
linear rates is crucial. Our proof idea breaks down for non-linear rates. The fact 
that non-empty closed intervals are used for specifying the transitions of A is not 
important. However, the fact that we have rectangular constraints is important. 

We now wish to argue that we can easily cope with networks of lazy hybrid 
automata in which the component automata communicate by synchronizing on 
common actions. 

Let V be a finite set of agent names with u,v ranging over V. We define 
a product lazy hybrid automaton to be a structure A-p = Y\ ue p A u where 
A u = (Q u , Act u , qf n , , D, {Pq} q eQ u , B, — > u ) for each u in V. For convenience, 

we will write T S u instead of T S. a u to denote the transition system over the reach- 
able configurations of A u as defined in section 2. The operational behavior of 
n uev A u is given by the transition system denoted as TS-p and is defined in the 
obvious way; it is just the usual synchronized product of the transition systems 
{TS'u}. The only twist is, in line with our discrete time semantics, all the compo- 
nents must move during a transition. If it is an a-move, then all the components 
that have a £ Act u must make an a-move while the remaining components must 
make a r-move. In a r-move, all the compoents must make a r-move. 

It is a routine exercise to establish that C st {Autp) and C act (Ap) can be 
computed effectively. 

5 Conclusion 

We have formulated here the class of lazy hybrid rectangular automata. These 
are basically linear rectangular hybrid automata but where each automaton is 
accompanied by the delay parameters {g, S g , h, Sh}- Our main result is that the 
discrete time behavior of these automata can be effectively computed if the 
allowed ranges of values for the variables are bounded. We have not outlined the 
verification problems for lazy rectangular hybrid automata that can be settled 
effectively. It should be clear however that we can model-check the discrete time 
behavior of our automata against a variety of linear time and branching time 
temporal logic specifications. 

We believe that associating non-zero bounded delays with the sensors and 
actuators is a natural assumption and it cuts down the the expressive power 
of hybrid automata. We also feel that it is useful to focus on the discrete time 
behavior of hybrid automata. Finally, there is some hope that larger classes of 
lazy hybrid automata may turn out to be tractable in terms of their discrete 
time behaviors. A related intersting problem is to determine the border between 
the decidable and undecidable in the context of lazy hybrid automata. 




Lazy Rectangular Hybrid Automata 



15 



References 

[1] R. Alur and D.L. Dill. A theory of timed automata. Theoretical Comp. Sci., 
126(2) : 183— 235, 1994. 

[2] R. Alur, T.A. Henzinger, G. Lafferriere, and G.J. Pappas. Discrete abstractions 
of hybrid systems. Proc. of the IEEE, 88:971 984, 2000. 

[3] E. Asarin, O. Bournez, T. Dang, and O. Maler. Reachability analysis of piecewise- 
linear dynamical systems. In Hybrid Systems: Comp, and Control, LNCS 1790, 
pages 20-31. Springer- Verlag, 2000. 

[4] V. Gupta, T.A. Henzinger, and R. Jagadeesan. Robust timed automata. In HART 
97, LNCS 1201 , pages 331-345. Springer- Verlag, 1997. 

[5] T.A. Henzinger. Hybrid automata with finite bisimulations. In 22nd ICALP, 
LNCS 944 1 pages 324-335. Springer- Verlag, 1995. 

[6] T.A. Henzinger. The theory of hybrid automata. In Proc. of the 11th Ann. Symp. 
on Logic in Comp. Sci., pages 278-292. IEEE Comp. Society Press, 1996. 

[7] T.A. Henzinger and P.W. Kopke. Discrete-time control for rectangular hybrid 
automata. Theoretical Comp. Sci., 221:369-392, 1999. 

[8] T.A. Henzinger, P.W. Kopke, A. Puri, and P. Varaiya. What’s decidable about 
hybrid automata? J. of Comp, and Sys. Sci., 57:94-124, 1998. 

[9] T.A. Henzinger and J.-F. Raskin. Robust undecidability of timed and hybrid 
systems. In HSCC 00, LNCS 1790, pages 145-159. Springer- Verlag, 2000. 

[10] Y. Kesten, A. Pnueli, J. Sifakis, and S. Yovine. Integration graphs: A class of 
decidable hybrid systems. In R.L. Grossman, A. Nerode, A. Ravn, and H. Rischel, 
editors, Hybrid Systems, LNCS 736, pages 179-208. Springer- Verlag, 1993. 

[11] J. McManis and P. Varaiya. Suspension automata: A decidable class of hybrid 
automata. In 6th CAV, LNCS 818, pages 105-117. Springer- Verlag, 1994. 

[12] J. Ouaknine and J. Worrell. Revisiting digitization, robustness and decidability 
for timed automata. In 25th LICS, pages 198-207. IEEE Press, 2003. 

[13] A. Puri and P. Varaiya. Decidability of hybrid systems with rectangular differen- 
tial inclusions. In 6th CAV, LNCS 818, pages 95-104. Springer- Verlag, 1994. 

[14] A. Tarski. A Decision Method for Elementary Algebra and Geometry. University 
of California Press, 1951. 




Affine Hybrid Systems 



Aaron D. Ames and Shankar Sastry 

University of California, Berkeley CA 94720, USA 
{adames , sastry}@eecs .berkeley . edu 



Abstract. Affine hybrid systems are hybrid systems in which the dis- 
crete domains are affine sets and the transition maps between discrete 
domains are affine transformations. The simple structure of these sys- 
tems results in interesting geometric properties; one of these is the notion 
of spatial equivalence. In this paper, a formal framework for describing 
affine hybrid systems is introduced. As an application, it is proven that 
every compact hybrid system H is spatially equivalent to a hybrid sys- 
tem Hid in which all the transition maps are the identity. An explicit 
and computable construction for Hid is given. 



1 Introduction 

This paper introduces affine hybrid systems. Affine hybrid systems are hybrid 
systems where the discrete domains are affine sets, and the transition maps be- 
tween discrete domains are affine transformations. This definition differs from 
other definitions of hybrid systems that have been proposed [9] , but the under- 
lying ideas involved in the definition of affine hybrid systems have been seen in 
the literature [6,7]. We give a formal framework to these ideas. 

Affine hybrid systems are simple, and it is this simplicity that allows us to say 
some useful things about them. The structure of affine hybrid systems contains 
a wealth of intrinsic information. Affine sets can be described in terms of matrix 
inequalities, and affine transformations are characterized by elements of SE(n). 
In this paper, we use the geometric information intrinsic in affine hybrid systems 
to develop the idea of spatial equivalence between an affine hybrid system H and 
an affine hybrid system G. 

In the literature on hybrid systems, it typically is assumed that all of the 
transition maps of a hybrid system are the identity; all switched systems are 
essentially hybrid systems where the transition maps are the identity [4,10]. 
This assumption is very restrictive; some of the simplest hybrid systems do not 
satisfy this assumption, e.g., the hybrid system T 2 constructed in Example 2.1 of 
this paper. For this reason, it is desirable to find a way to bridge the gap between 
hybrid systems where all the transition maps are the identity and hybrid systems 
where this is not the case. 

Given an affine hybrid system H, we would like to construct an affine hybrid 
system H;d such that all of the transition maps are the identity. We also would 
like this affine hybrid system Hid to be as similar to H as possible. In what way 
should these two affine hybrid systems be considered similar? Spatial equivalence 
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is introduced as a way to consider an affine hybrid system H as similar to an 
affine hybrid system G. Spatial equivalence can be thought of in an intuitive 
manner (see Figure 4 for a visual interpretation). Replace each edge of H by 
a sequence of edges and domains with vector fields such that if we “start” at 
the source of the edge, the target of the edge will be reached in some time. If 
the affine hybrid system obtained by appending these edges, domains and vector 
fields to H is G, then H is spatially equivalent to G. A formal definition of 
spatial equivalence will be given in Section 5, but having this intuitive picture 
in mind will be helpful. 

An affine hybrid system H is compact if each of its domains is compact. The 
main theorem of this paper is: 

Main Theorem. Every compact affine hybrid system H is spatially equivalent 
to an affine hybrid system Hid in which every transition map is the identity. 
Moreover, H;d is computable. 

This paper begins by introducing, in Section 2, the definition of an affine 
hybrid system. Sections 3 and 4 present some results in affine geometry that are 
necessary for the proof of the Main Theorem. In Section 3, given two (n — 1)- 
dimensional affine sets X and Y = RX + p, for ( R,p ) £ SE(n), we determine 
conditions on X, R and p such that there exists an n-dimensional affine set 
with X and Y as faces. When these conditions are satisfied, we find a closed 
form solution for a set S which has A' and Y as faces. This closed form solution 
allows us later to compute Hid- When there does not exist an affine set S with 
X and Y as faces, an admissible sequence of faces, X = Zq, Z-[, .... Zf, = Y, is 
introduced; it is used to construct a sequence of affine sets, Si, S ?, ..., Sh, where 
X is a face of Si and Y is a face of Sk- Admissible sequences of faces are used in 
Section 4 to generalize the results of Section 3 by showing that if Y = RX + p, 
for (R,p) £ SE(n), there is an admissible sequence of faces 

Zo = -A ■Zlj *•*! ^^n(n-l)+l) Z ™ n {n-l)+2 = Y. 

In Section 5, the results in affine geometry that were introduced in Sections 3 
and 4 are used to prove the Main Theorem and give an explicit construction for 
Hi d . This is done by using the admissible sequence of faces found in Section 4 
to define the domains of the hybrid system Hid- 

2 Affine Hybrid Systems 

This section introduces the notion of an affine hybrid system. An affine hybrid 
system consists of the following data: a set of discrete states, domains, edges 
and vector fields. The discrete states provide a way to index the domains. The 
domains are affine sets, i.e., sets that are affinely constrained. The edges provide 
a relationship between two faces of two domains; each edge has a source which 
is the face of a domain and a target which is also the face of a domain. It is 
required that there exists an affine transformation between the source and the 
target of each edge; thus each edge gives rise to a transition map, which is an 
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affine transformation from the source of the edge to the target of the edge. The 
set of vector fields is a collection of vector fields that are Lipschitz on each 
domain. 

Before we formally introduce the definition of an affine hybrid system, we 
will describe each of the components of its definition. This section is concluded 
by solidifying the concepts introduced through an example: the torus T 2 . This 
example also will be important later in the paper. 

2.1 (Discrete states). Let Q C Z denote the set of discrete states. This set is 
finite, and the number of discrete states is given by \Q\. For simplicity, typically 
we let Q = {1, ..., m}. 

2.2 (Domains). The set of domains is the set D = {Di}i & Q, where each Di C R” 

is an n-dimensional affine set, i.e., a set that is affinely constrained. For each set 
Di, there exists a matrix A,; £ and a vector a-, £ K fci such that 

x £ Di <t=> AiX + a* > 0. 

We say that Di is determined by the affine constraints AiX + at . 

The boundary of Di can be written as the union of hi affine sets of dimension 
n — 1 called the faces of Di. The faces of Di can be indexed by introducing the 
indexing set, 

Fi = {l,...,ki}, i £ Q. 

The j th face of Di is denoted by Facej(Dj), where j £ Fj. We can pick an 
indexing of the faces of Di by letting Facej (£>.;) be the affine set determined by 
the j th row of Aj. More precisely, if (Aj)j* is the j th row of A, and (m)j is the 
j th entry of a^, then 

x £ Face., (A) ^ (-(A*),-*) iE+ (-(a;),) - °‘ (1) 

We can define 

Aij = (-(to,*)’ ° ij= (-K)J’ 

so x £ Facej(D,) if and only if A,jX+a,y > 0. Therefore, Facej(Dj) is determined 
by the affine constraints A*jX + a^. 

2.3. For a set U with U = n"=i denote the projections on each of the factors 
of U by TTi : U — > Ui . 



2.4 (Edges). Define the set of edges as a set 

EC {((*, j), (M))}(i,j)eQ X Q, (k,i) eFixFj, 

satisfying the condition that for each e £ E, there exists a map T e {x) = R e x+p e , 
with ( R e ,p e ) £ SE(n), such that 

T e (Face,^ ( e ) (D Vl ( e ) ) ) Face-g^g^-D^^g)). 



(2) 
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In other words, an edge defines a relationship between two faces of two affine 
sets and an affine transformation between these faces. 

More concretely, an element e £ E then has the form 

e = ((bj)>(M)), {i,j)€QxQ, (k,l) £ FiX Fj, 

so 7Ti(e) = i, 7T2(e) = j, n^e) = k and 7r4(e) = l. Condition (2) allows us to 
write 

T e (Face fc (A)) = -R e Face fc (A) +Pe = Face ; (A)- 



2.5. Given an edge e £ E, the affine transformation T e {x) = R e x + p e from 
Face 7r3 ( e )(D 7r i( e )) to Fac e n4 r e AD n2 / e \) is called the transition map. The set of 
transition maps is the set T = {T e } e& E- 



2.6 (Vector field). A set of vector fields is a set V = {V)}j £ Q where V) is a 
Lipsclritz vector field when restricted to the domain A. The flow of V) on A is 
denoted by <pi(t,x) for x € Di. 



Definition 2.1. An affine hybrid system is a tuple H = (■ Q,D,E , V). 



Note 2.1. From this point on, for the sake of brevity, we will refer to “affine 
hybrid systems” as “hybrid systems”. 



2.7. If for some e £ E, T e (x ) = x, then we say that the transition map associated 
to the edge e is the identity map. This implies that 

Face 7r3 ( e )(.D 7 r 1 ( e )) Fac e 7r4 ^ e ^(D^. 2 ^ e ^). 

A very special class of hybrid systems are hybrid systems in which every transi- 
tion map is the identity, and we denote such hybrid systems as Hij. 



Example 2. 1 (The torus: T 2 ). We will construct a hybrid system called the torus, 
which we will denote by T 2 (see Figure 1). The torus is given by one discrete 
state Q j2 = {1}. The domain Dj 2 = {(x,y) : x £ [0,1], y £ [0,1]} is given by 
the affine constraints 



Aj 



\ , T 2 


( 1 o \ 
-1 0 


( x\ 


(°\ 

l 


) +4 = 


0 -1 


U) + 


l 




l 0 1 ) 
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Applying (1), the affine constraints for Facei (Df 2 ), Face 2 (Dj 2 ), Faces {D J”) and 
Face 4 (Dj ) are given, respectively, by 




£ t2 consists of two edges: e\ = ((1, 1), (2, 1)) and e 2 = ((1, 1), (3,4). In other 
words, ei is a relation between the top and bottom of the square and e 2 is a 
relation between the right and left side of the square. The associated transition 
maps are 



T ei ( X ,y) = (*) + ( q 1 ) , T e2 (x,y) = + (_°j) • 

Finally, V^ 2 is any vector field, Lipshitz on Df 2 . 

The advantage of defining the edges as a relationship between the faces rather 
than a relationship between the domains can be seen in this example. Although 
the expression for the edges is more complicated, the end result is a simpler 
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definition of the hybrid system; in other references, the torus is defined with two 
discrete states [9]. 



3 Affine Sets 



Given two (n — l)-dimensional affine sets A and Y = RX+p, for (R,p) £ SE(n ), 
is it possible to find an affine set S with X and Y as faces? Clearly the answer 
to this question is no for an arbitrary element of SE(n), but it is yes if X is in 
the “proper position” and R and p satisfy certain assumptions. The purpose of 
this section is to find a closed form solution for the affine constraints defining a 
set S with X and Y as faces, when the assumptions on X, R and p are satisfied. 
This result is important because it makes the later propositions and theorem of 
this paper computable by way of this closed form solution. We also will use this 
formula repeatedly in order to compute examples, beginning with an example 
at the end of this section. For more detailed proofs of the results presented in 
this section, see [1,2]. 



3.1. First, recall some important facts and terminology regarding affinely con- 
strained sets. We define a face of an n-dimensional affine set X, denoted by 
Facei(A) for i = 1, ...,k (where k is the number of faces), as a subset dX such 
that there exists a lryperplane Hi where Hi (~l dX = Face, (A). This hyperplane 
is called the lryperplane defining Facei(A). If X is determined by the affine con- 
straints Ax + a , and if we define the Face j (A) as the set determined by the affine 
constraints 



AiX + a* 




then the defining lryperplane, Hi, is given by Hi = 1 aijXj +ai = 0}. If the 

smallest number of affine constraints that determine A is k, then A has k faces. 
Note that in this case it is always possible to define A in terms of more that k 
affine constraints, but never less. 



Proposition 3.1. Let X be an affine set of dimension n— 1 in R™, and assume 
that X C { Xi = 0}. IfY = X + p, with Pi ^ 0, then there exists an affine set S 
such that X andY are both faces of S. Moreover, there is a closed form solution 
for the affine constraints that determine S. 

Proof. If A C {a;,; = 0} is an (n — l)-dimensional affine set with k faces, then 
the affine constraints defining A can be put in the form 



/an • • 


• ai,i_i 0 


ai,i+i • • 


’ O'ln \ 
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Since Pif^ 0, if we define 



1 

Cfc — 

Pi 



yi Pj a kj > 

i = i 
j / * 



the affine constraints for the set S are given by 



/an • • 


• ai,i_i 


Cl 


ai,;+i • • 


’ O'ln ^ 




( X\\ 




( ai N 


aki ■ ■ 

0 •• 
V 0 •• 


' ak,i— 1 

• 0 
• 0 


Cfc 

sign (pi) 

-sign(pi) 


ak,i+i ' ■ 

0 •• 
0 •• 


CLkn 

• 0 
• 0 ) 




\X n J 


+ 


0 

\sign {pi)pij 



It can be verified easily that X is the face of S given by intersecting S with the 
lryperplane {a;,: = 0}. Similarly, Y is the face of S given by intersecting S with 
the hyperplane {xi — pi = 0}. □ 



3.2. Throughout this paper, we will use angle to refer to a scaler with values in 
[— 7T,7r). For n > 2, Givens rotations (see [5,8]) are n x n matrices of the form 



column i column j 



i 



i 



cos 6 



Pij(0) = 



— sin 6 



sin 6 



cos 6 



i / 



Givens rotations are important because, for every R £ SO{n) with n > 2, 
there exists n(n — l)/2 angles % £ [— 7r,7r) such that R = n"=i llj=i+i 
(cf. [3]). Moreover, there is a closed form solution for . Therefore, understand- 
ing the effect of applying an element of SO{n) to an affine set is equivalent to 
understanding the effect of applying a Givens rotation. 



Proposition 3.2. Let X be an affine set of dimension n— 1 in K™, and assume 
that X C {xi = 0} fl {xj > 0}. If Y = Pij(6)X , with 6 £ (0,n), then there 
exists an affine set S such that X and Y are both faces of S. Moreover, there is 
a closed form solution for the affine constraints that determine S. 
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Fig. 2. Left: The sets X and Y in Example 3.1. Right: The set S with X and Y as 
faces. 

Proof. If X C {xi = 0} fl {xj > 0} is an (n — l)-dimensional affine set with k 
faces, the affine constraints defining X can be written as 



f an • 


• ai,i_i 


0 ai,i+i • 


O'ljj 1 Qrl j 


aij+i • 


■ ai n \ 




( x l\ 
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The affine constraints for the set S with X and Y as faces are given by 



column i 

i 



column j 

i 



/an • ■ 


• ai,i_i 
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0 •• 
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V o ) 



( 6 ) 



• ( 7 ) 



It can be verified easily that X is the face of S given by intersecting S with the 
lryperplane {a:, = 0}. Similarly, Y is the face of S given by intersecting S with 
the lryperplane {cos to* + sin Oxj = 0}. □ 



Example 3.1. Consider the set X = {(a;, y) : x = 0, y G [0, 1]} and Y = Pi 2 (f )X. 
Since X C {a; = 0} fl {y > 0}, we can apply Proposition 3.2 to determine an 
affine set S with X and Y as faces. The affine constraints for X are given by 
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where these affine constraints are in the same form as (6). Applying (7) gives 
the affine constraints for S as 



/ (1 - V2) 1 

-( 1 -^ 2 ) -1 



l 



\ , N 


(°\ 




i 


\y) + 


0 


/ 


W 



Or this set is given by the constraints that y < —(1 — \fT)x + 1, y > — x and 
x < 0. The remaining constraint, that y > —(1 — \/2)x , is satisfied when the 
other three constraints are satisfied. The set S is exactly the set that we would 
have hoped for (see Figure 2). 



4 Admissible Sequences 

For two (n — l)-dimensional affine sets X and Y = RX + p, in general it is 
not true that there exists an n-dimensional affine set with X and Y as faces. 
When there is not an affine set with X and Y as faces, the question is: does 
there exist a sequence of affine sets where the first affine set has A as a face, the 
last affine set has Y as a face, and any two adjacent affine sets in the sequence 
share a common face? In this section, it will be shown that for any set X and 
Y = RX + p , there exists a sequence of affine sets of this form; these sequences 
will be essential to the proof of the Main Theorem. The results of the previous 
section allow each of the affine sets in the sequence to be computed. Detailed 
proofs of the results of this section can be found in [1,2]. 

Definition 4.1. Two (n — 1)- dimensional affine sets, X and Y, are admissible 
faces if there exists an n-dimensional affine set S(A, Y) with X and Y as faces. 



4.1. If E(X, Y) is an affine set, for ( R,p ) £ SE(n), there are the following 
properties 



E(A,Y) = E(Y,A), 

RS(X, Y)+p= E (RX + p, RY + p ) , 

where E (RX +p,RY + p) is an affine set with RX +p and RY +p as faces. 



4.2. We have shown in Proposition 3.1 that if X C {xi = 0} and Y = X + p 
with pi 0, then X and Y are admissible faces; we can take E(A, Y) to be the 
affine set given by the affine constraints in (5). Similarly, by Proposition 3.2, if 
X C {xi = 0} (~l {xj > 0} and Y = Pij(0)X, for 0 € (0,7r), then X and Y are 
admissible faces; we can take E(A, Y) to be the affine set given by the affine 
constraints in (7). 
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Fig. 3. Left: the sets X, Z\ and Y in Example 4.1. Right: the affine sets Si = E(A', Zf) 
and Si = S{Zi,Y). 



Definition 4.2. A sequence Z 0 , Z x , Z^ of {n — 1)- dimensional affine sets is 
an admissible sequence of faces if there exists affine sets, 

^(Zq,Zi), S(Zi,Z 2 ), ..., 3(Z fc _i, Z fc ). 



Proposition 4.1. Let X be an (n — 1) -dimensional affine set and Y = X + p. 
Then there exists an admissible sequence of faces Zq = X, Z\, Z 2 = Y . 

Proposition 4.2. Let X be a compact (n— 1) -dimensional affine set with n > 3, 
and Y = Pij(9)X. Then there exists an admissible sequence of faces 

Z 0 = X, Zi , ..., Z 9 , Zh = Y. 



Remark fl. In the case where n = 2, an obvious modification of Proposition 4.2 
gives an admissible sequence of faces Z 0 = X , Z x , ..., i? 4 , Z 5 = Y. Throughout 
the rest of the paper, we will assume that n > 3. All of the results are applicable 
to the case where n = 2, with the obvious modifications. 



Theorem 4.1. Let X be a compact (n — 1)- dimensional affine set, and Y = 
RX +p, with ( R,p ) £ SE(n). Then there exists an admissible sequence of faces 

Zo = X, Z \, ..., Zu n ( n _ i) +1 , Z±i n ( n _ 1 ) + 2 = 1 . 



Example f.l. Let X = {(x,y) : x £ [0,1], y = 0}, p = (2,0), and Y = X + p, 
i.e., Y = {(x,y) : x £ [2,3], y = 0}. Therefore, X and Y are given by the affine 
constraints: 
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respectively. It is clear that there is not a single n-dimensional affine set with X 
and Y as faces. This is evident in the fact that the assumptions of Proposition 
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3.1 are not satisfied; X C {y = 0}, but p 2 = 0. By Proposition 4.1, we can find 
a sequence of admissible faces Z 0 = X, Z \ , Z 2 = Y, and the corresponding affine 
sets Si = S (X,Zi) and S 2 = S (Zi,Y). 

Define a = (2, 1) and b = (0, —1), then a + b = p and Z\ = X + a = {(&, y) : 
x £ [2,3], y = 1}. We can let Si = S(X, Z{), which is given by the affine 
constraints in (5). Since Z\ = Y — b, S 2 = a(Zi,Y) = 3 (Y — b,Y) = 3 (Y,Y — b). 
Because Y C {y = 0} and b- 2 ^ 0, applying Proposition 3.1 gives the affine 
constraints for S 2 . Therefore, applying (5) to the affine constraints defining X 
and Y gives the affine constraints defining S± and S 2 as: 



( 1 
-1 2 


( x \ 


1 




/ 1 0 \ 
-1 0 


( x \ 


(~ 2 \ 

3 


0 1 


u + 
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0 1 


u + 


0 


V 0 -v 




v) 




l 0 -l) 




^ 1 / 



respectively. For a visual interpretation of these results, see Figure 3. 



5 Spatial Equivalence 

In this section we use the sequence of admissible faces found in Section 4 to 
prove the Main Theorem of this paper. In the proof of this theorem, the hybrid 
system H;<j is constructed. This allows H;d to be explicitly computed, as will be 
seen through an example following the proof of the Main Theorem. 

5.1. A hybrid system H = (Q u , D H , E H , V H ) is said to be spatially equivalent 
to a hybrid system G = (Q , D G , E G , V G ) if the following conditions hold: 

1. |Q g | > |<3 H | = m. 

2. For every i < to, Dj 1 = D G and for i > to, there exist admissible faces X t 
and Yi such that D G = 3 (JQ,Y;). 

3. For every edge e £ E H there exists a sequence of k discrete states 1), ..., 
y{k) > to and edges r/i , ...,Tjk+ i £ E G such that 



^»7i(Face 7r 3 (e)(-Z5 7ri (e))) — 



L ^(i) ’ 



T V2 (Y v (i)) — 



’^'vk+i(Yv{k)) Fac 

In the special case where k = 0, we require pi to be an edge such that 

T Vl (Face^g (e) (D^ (e) ) ) = Fac e„ 4(e) (I?“ (e) ). 

4. For i < to, E H = V G , and for i > m, V G is a vector field such that 
<j> G (l,Xi) = Yi, where 4> G (t,x) is the solution to V G . 
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Remark 5.1. Note that spatial equivalence is not an equivalence relation. The 
term “equivalence” is used in order to stress the equivalence of the qualitative 
behavior of H and G, when H is spatially equivalent to G. Although it might be 
appropriate to replace “spatial equivalence” by a term such as “spatial embed- 
ding” , the authors are concerned that this term would not stress the behavioral 
similarities of the two hybrid systems. 



Definition 5.1. A hybrid system H is compact if each of its domains is com- 
pact. 



Theorem 5.1. If H is a compact hybrid system, then H is spatially equivalent 
to a hybrid system H; d such that every transition map is the identity, i.e., T v = id 
for every rj £ E Uid . Moreover, Hy is computable. 

Proof. In order to prove this theorem, we will explicitly construct the hybrid 
system H; d . First, we define 

^:={ee£ H : T e # id}, 

which is the set of edges such that the associated transition map is not the 
identity. If |-E?f| = k, then we can write the elements of £yf as ei,...,efc (by 
arbitrarily indexing them). For simplicity of notation, define the functions 

/( n ) = y n ( n - !) + 2 > 
g(in , n, i) = (i — 1 )f(n) + m + 1, 

which will be used throughout the course of the construction. 

Construction of H id 

Q Hid : If Q h = {1, then define Q Hid = {1, ..., in + kf(n)}, with k = |-Eyf|. 

D Uid : For i < m, define £> ! Hid = Df 1 . If AiX + a.j are the affine constraints 
determining Dj 1 , then we also let AiX + ai (with the order of the rows main- 
tained) be the affine constraints determining D^' d . In particular, this implies 
that Face : ,(D* Iid ) = Facej(I>P). 

Now we can construct the other domains of Hi d . For every edge ei £ E™, 
1 < i < k, the transition map is given by T ei (x) = R ei x + p ei , or we have 

Face 7r4 ( ei )(D 7r2 ^ ei j) = i? ei Face 7r3 ( ei )(D 7ri ^ ei 'j) + p ei . 

Now by Theorem 4.1 we have the following admissible sequence of faces 
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Setting 



X, 



g(m,n,i ) 



= Z, 



0 ’ 



Y yi ~y 

y *-g(m,n,i)-\- 1 1 g(m,n,i) ? 



Xg(m,n,i)+f(n)- 1 ^fl(m,n,i)+/(n)- 2, 



Zf(n) ~ Yg(.™,n,i)+f(n)- 1, 



D , .■ 5 ^ r g(m,n,i)+j') , 1 ^ 7 ^ A:, Q < j — f ( 77 ) 1 - 



define the domains 

}H id 

' g(m,n,i)+j 



It can be verified that for these values of i and j, g(m , n, i) + j takes all values 
from to + 1 to m + kf(n), inclusive, and with no repeats. 



E Uid : If e £ E u and e ^ E™, then the associated transition map is T e = id. 
So define an edge 77(e) £ E Uid to be 77(e) = e. It follows that T v ( e ) = id. 
If e £ E™, then e = e, for i £ {l,...,fc}. We can now define a set of edges 
/7i(ei), 772(6*), ...,?7/(n)+i(ei) £ E Hid as follows: if we index the faces of D^^ ni)+j 
such that 

Xg(m,n,i)+j — Facei (-D g ( m „ j) +J )) }g(m,n,i)+j = FaC e 2 (^ g (m,n,i)+j^ 

then we define 

m (e*) = ((TTi(ei), 3(777, ?r, i)), ( 7 r 3 (e,:),l), 

772 (e*) = {(g(m,n,i), 3(777, n,i) + 1), (2,1)), 



77 j ( ei ) = ((3(771,77,*) + j - 2 , 3(777., 77 , i ) + j - 1 ), ( 2 , 1 )), 

Vf ( n )( e i ) = (( 3 ("i, 77 , 7 ) + f (77) - 2,3(777,77,7) + / (ti) - 1 ), ( 2 , 1 )), 
? 7 /(n)+i(e*) = ((3(777,77,7) + /(n) - l, 7 r 2 (e i )),( 2 , 7 r 4 (e i )). 

The associated transition maps are 

-^i(ei) • Face^ ( e ^) (-^ 7 ri (e i )) “ ^ ■Zg{m,n,i)i 

r ^' r )j{ e i) * Fp(m,n,i)+j— 2 ^ -^p(m,n,i)+j — 1 5 1 ^ j — 

^77/(„) +1 (e<) • — >■ Face^ ( ei ) ( D ^ 2 ( e . ) ) , 

so 

Face 7 r 3 ( e .)(D 7ri ( e .)) ■ Zg { m , n , i )'> 

Fp(m,n,i)+j— 2 — 1 5 ^ ^ J — /(^)? 

F^( m ,n,i)+/(n) — 1 = Face 7r4 ( ei ) (-^ 7r2 ( e< )) • 
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By definition, it follows that 

^rn (e») = ^»)2(ei) — • • • — T. ij /( „) + i(ei ) = id. 

If we apply this construction to every edge in E u the result is the set E Hid . It 
is clear that for every edge in 77 £ E Uid , T v = id. It also follows that |E Hid | = 
|£ H | + 2/(n). 



U Hld : If i < m, define V' Hid = V) H . If « > m, then where X,; 

and Yi differ by an element of SE(n), i.e. , Yi = QiXi + qi. Using this, we define 

4>f id (t, x ) = (1 - t)x + t{QiX + qi), V™- ld (x) = 
and we have the property that (jy* id {l,Xi) = QiXi + qi = Yi. 

To conclude the proof we note that in the process of constructing Hid we 
have shown that H and H;d satisfy properties 1-4 of Paragraph 5.1, hence H is 
spatially equivalent to H;d- □ 

Remark 5.2. Note that the hybrid system we constructed in the proof of The- 
orem 5.1 is not unique. Moreover, the number of discrete states given in the 
construction is not necessarily the smallest number of discrete states needed to 
construct a spatially equivalent hybrid system. For example, if for every edge 
e £ E h , the faces Face^^D^^) and Fac e W 4 ( e )(T ) ^( e )) are admissible (see 
Definition 4.1), then we can construct a hybrid system Hid, spatially equivalent 
to H, with Q Hid = {1, ...,m + k}. This will be the case in the following example. 



Example 5.1. We will construct T 2 d , or a hybrid system spatially equivalent to 
T 2 (see Example 2.1) where every transition map is the identity. In this case we 
have two edges, ei, e 2 £ E T , and | E^ \ = 2. First, note that we can define T 2 d in 
terms of fewer than the number of discrete states given in the proof of Theorem 

5.1 because Fac e 7 r 3 (e i )(F ) ^ I l(e .)) and Fac e, r 4 (e i )(F ) ^(e,)) are admissible faces for 
*=1,2 (see Remark 5.2). 

Set Q T ‘ d = {1,2,3}, and let Dj id = Dj , which we defined in Example 2.1. 
To construct D 2 id and D 3 id , note that 

Face 2 (Hf) = Facei(T)f ) + (j) . Face 3 (uf ) = Face 4 (Df) + ■ 

Since Facei(Z?} 2 ) C {x = 0} and Face^-D]" 2 ) c {y = 0}, by applying Proposition 

3.1 to the affine constraints Aj t x + a 41 and A\ a x + o} 4 , given in equations (3) 
and (4), it can be verified that 

H(Facei(Df ),Face 2 (£>f)) = H(Face 4 (Df), Face 3 (Df)) = D? . 
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Fig. 4. Left: T 2 . Right: lf d . 

Now as in the proof of Theorem 5.1, define 

X 2 = Face 2 (£>f ) , F 2 = Facer (£>f ) , 

X 3 = Faces (D ? ) , F 3 = Face 4 (D ? ) , 

then 

£ 2 fd = S(X 2 ,y 2 ) = Df, Df« = H(X 3 ,y 3 ) = Dj\ 

and D T 

Now we will determine the edges of T? d . As in the proof of Theorem 5.1, 

T 2 T 2 

index the faces of D 2 id and D 3 id such that 

X 2 = Facer (£> 2 ‘ d ) , Y 2 = Face 2 (Z? 2 ‘ d ) , 

X 3 = Facer (B? ) , Y 3 = Face, (H 3 ' d ) , 

and for the two edges er, e 2 £ E t2 , define 

Vi(ei) = ((1,2), (2,1)), ?7r(e 2 ) = ((1, 3), (4, 1)), 

? ?2 (er) = ((2, 1), (2, 1)), , ?2 (e 2 ) = ((3, 1), (2, 3)). 

Set E Tid = {??r(er), ?? 2 (er), ??r(e 2 ), ? 72 (e 2 )}, and note that the corresponding tran- 
sition maps T m(ei) = T^ 2 ( ei ) = r m(e2 ) = T^ 2(e2 ) = id (see Figure 4). 

Finally, define V 1 Tid = F 1 T ’" and 

V^{x) = (" q 1 ), ^ 3 Tfd (^)=(_ 0 1 )- 

So, F T ‘ ?d = {F 1 Tld ,y 2 Tid ,y 3 Tid }. This completes the construction of T? d . 
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Abstract. In this paper we present an abstraction method for nonlin- 
ear continuous systems. The main idea of our method is to project out 
some continuous variables, say z, and treat them in the dynamics of the 
remaining variables x as uncertain input. Therefore, the dynamics of x is 
then described by a differential inclusion. In addition, in order to avoid 
excessively conservative abstractions, the domains of the projected vari- 
ables are divided into smaller regions corresponding to different differen- 
tial inclusions. The final result of our abstraction procedure is a hybrid 
system of lower dimension with some important properties that guaran- 
tee convergence results. The applicability of this abstraction approach 
depends on the ability to deal with differential inclusions. We then focus 
on uncertain bilinear systems, a simple yet useful class of nonlinear dif- 
ferential inclusions, and develop a reachability technique using optimal 
control. The combination of the abstraction method and the reachabil- 
ity analysis technique for bilinear systems allows to treat multi-affine 
systems, which is illustrated with a biological system. 



1 Introduction 

Recent developments in embedded control systems have motivated much research 
on automated verification of continuous and hybrid systems. For systems involv- 
ing non-trivial continuous dynamics (described by differential equations), exact 
and approximate reachability analysis methods have been developed [16,22,9,3, 
20,24]. Even though these methods have been used to treat interesting case stud- 
ies, the complexity of reachability computations currently limits the application 
to small size systems. In order to scale to larger systems, abstraction methods 
have been investigated (see [9,32,2,31,26]). Roughly speaking, abstraction is a 
general approach allowing to deduce properties of a system by analyzing a more 
abstract and, in general, smaller system (see [10,8] and references therein). Some 
of the abstraction methods for hybrid systems, inspired by techniques from pro- 
gram analysis, aim at extracting (exactly or approximately) a finite-state model 
from a continuous/hybrid system while the other exploit the structure of the 
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system in order to reduce it into a system of smaller dimension which preserves 
the properties of interest. In addition, the methods based on approximating con- 
tinuous systems by hybrid systems with simpler continuous dynamics [28,18,4], 
that we call ‘hybridization-based’ methods, can also be viewed as an abstraction 
approach. 

In this work, we propose an abstraction method for dimension reduction, which 
is along the lines of the hybridization-based approach. Our first observation is 
that in many practical systems, the properties to verify involve only a subset 
of variables, and the other variables may not need to be analyzed with great 
accuracy. The main idea of our method is to project out some continuous vari- 
ables, say 2 , and treat them in the dynamics of the remaining variables x as 
uncertain input. Therefore, the dynamics of x is then described by a differential 
inclusion. In addition, in order to avoid excessively conservative abstractions, 
the domains of the projected variables are divided into smaller regions corre- 
sponding to different differential inclusions. The final result of our abstraction 
procedure is a hybrid system of lower dimension with some important proper- 
ties that guarantee convergence results. However, this abstraction method does 
not solve the verification problem by itself. The success depends on the ability 
to deal with differential inclusions. We thus focus on the reachability problem 
for uncertain bilinear systems, a simple yet useful class of nonlinear differen- 
tial inclusions. The combination of the abstraction method and the reachability 
analysis method for bilinear systems allows to treat multi-affine systems, which 
can be found in numerous applications in engineering, biology and economics. 
The rest of the paper is organized as follows. In Section 2 we present our abstrac- 
tion method and the convergence results. Section 3 is devoted to a reachability 
analysis method for uncertain bilinear systems, which is motivated by the ap- 
plication to multi-affine systems. This reachability technique uses results from 
optimal control. Section 4 contains an example of a biological system illustrating 
the theoretical results of the paper. 

2 Abstraction by Projection 

2.1 Basic Idea 

We consider a continuous system 



(x = f(x,z) 
= g(x,z) 



(1) 



where x € X C R n , z € Z C R m . We assume that the state space of the system 
is compact and the functions /, g are Lipschitz continuous. Given a vector x we 
use the notation x t to denote the i th component of x. 

Suppose that we want to reduce the dimension of the system (1) from n + m to 
n by projecting out the variables z. As in qualitative simulation, the first step 
of the abstraction is to partition the domain of z into disjoint regions, and in 
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each region the dynamics of x is approximated by a differential inclusion which 
is obtained from f(x,z ) by letting z to take any value in the region. In other 
words, the effect of z in the dynamics of x is modeled as uncertain external 
input. Let us now formalize this idea. 

We suppose that the domain of variable Zi is an interval It = [z iy Zi], and the 
domain of z is thus a box B = X\ x Z 2 x . . . x X m . We partition each interval Xi into 
k disjoint intervals 1 of the form {I 1 = \zj,zj), Xf = [ z 2 , zf), . . . ,Zf = [zf, zf]} 
such that zj = z t , zf = 2 , and for all j £ {1, . . . , fe — 1} zf = zf +1 . Therefore, 
the box B is partitioned into k m boxes, and we denote by B 1 with i £ N m the box 
ZJ 1 x ... x Z^ 1 . In the following we shall approximate the (n + m)-dimensional 
continuous system (1), referred to as the original system, by a hybrid automaton 
with n continuous variables. 

Each box B 1 corresponds to a location loc 1 of the approximating hybrid au- 
tomaton where the dynamics of x is approximated by the following differential 
inclusion: x £ F 1 (x) = { f(x,z ) | 2 £ B 1 }. The transitions of this hybrid au- 
tomaton correspond to the reachability relation between the boxes B 1 of the 
original system (1), which can be abstracted as follows. Note that since (1) is 
continuous, only transitions between adjacent boxes need to be considered. For 
our further developments, we need to introduce some additional notations. The 
boxes B 1 and B 3 are called adjacent if |j t — i,| < 1 for all i and i / j. We denote 
by d( i, j) = {i | 1 < i < m A j* ^ i,;|} the set of indices at which the components 
of i and j differ. We use d(B\ B •*) to denote the common boundary of the boxes. 
Given two adjacent boxes B 1 and B 3 , the condition for the transition from B 1 to 
B 3 , denoted by B 1 — > B 3 , is: 

£ d(B\B 3 ) Vi £ d(i,j) : (j; - i i)gi{x,z) > 0 (2) 

where gt denotes the i th component of g. The above condition says that the 
transition from a box to one of its adjacent boxes is possible if there exists at 
least one point on the common boundary of the two boxes at which the derivative 
of 2 points into the arrival box. As an example, for two adjacent boxes B 1 and B 3 
such that j,; = ij + 1 and j_, = i ; for all other j ^ i, the condition for the transition 
B 1 — > B 3 is gffx^z 1 ) > 0. Similarly, the condition for the transition B 3 — > B 1 is 
gi(x,z l ) < 0. Obviously, the condition (2) is not sufficient since it only implies 
that there exists a trajectory of the original system that goes from one box to 
an adjacent box. Hence, the resulting hybrid system is an over-approximation of 
the original continuous system. The approximating system where 2 is scalar is 
shown in Figure 1. 

2.2 Remedy Discontinuities 

It should be noted that the way we project out the variables 2 introduces dis- 
continuities in the derivative of the remaining variables x. For the sake of well- 
posedness, as in the sliding mode 2 approaches, we shall “convexity” the dynamics 

1 For simplicity of notation, we choose the same number of intervals for all z t . 

2 The literature of sliding mode control is vast, see for example [34,1]. 
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g(x,z 1 )>0 g(x,z 2 )> 0 g(x,z k > 0 




g(x,z 1 )<0 g(x,z 2 )< 0 g(x,z k x ) < 0 



Fig. 1. Hybrid automaton obtained after projecting out the variables z. 



of the abstract system in order to guarantee the existence of solutions [14] as 
well as an error bound of the approximation of the solution set. 

Let B 1 and B j be two adjacent boxes. Between the locations loc 1 and loc?, we 
add a location loc u whose continuous dynamics is defined as: x € F 1} (x) = 
colF 1 (x) , F* (x)} where co denotes the closed convex hull. The resulting system 
is illustrated in Figure 2. For brevity, we denote the approximating system of 
Figure 2 by 

x € T(x). (3) 




~^g(x, z { ) > o 



x £ F l (x, y) 




g(x,z x ) > 0 



g{x , z J ) >0 



x € F'd x, y ) 

9 (x,z i ) = 0 

g{x,z i ) < 0 " ^ g(x,z l ) < 0 




Fig. 2. Approximating hybrid system with upper semi-continuous dynamics. 



The above “convexification” provides the system (3) with an important property 
stated in the following lemma. 

Lemma 1. The multifunction T in (3) is one-sided Lipschitz and upper semi- 
continuous. 

The concepts of one-sided Lipschitz and upper semi-continuity are recalled in 
the proof of the lemma, presented in Appendix. 

To quantify the error between the original and approximating systems, we first 
define the size of the discretization of 2 . The diameter of the box B 1 is diam(B 1 ) = 
max {|2 — z'\ : z € B 1 A z' £ B 1 } where | - 1 is the Euclidean norm. Then, the size 
of the discretization of z is 5 Z = maxi{dza?n(S 1 )}. 

Lemma 2. Let 6 Z be the size of the discretization of z and Lf be the Lipschitz 
constant of the function f in (1). Let be a solution of (1). Then, for 

all t > 0 dH(x(t),F(x(t)) < LfS z where dn is the Haussdorf distance. 
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The above lemma gives a bound on the distance between the derivatives of the 
original system and the approximating system. Using the one-sided Lipschitz and 
upper semi-continuity properties of the approximating system, we can establish 
the following convergence result. 

Theorem 1. Let #(•) and z(-) be absolutely continuous functions satisfying (1). 
Then, there exists a solution x{-) of (3) such that for all t > 0 

\x(t) - x(t ) I < KO) - m\e Ct + ^(e ct - 1) 

where C is the Lipschitz constant of T , and A = Lf8 z is the bound on the 
distance between the derivatives estimated in Lemma 2. 

This theorem is proved in Appendix. The proof uses the assumption that T{x) 
is nonempty and takes closed convex values. 

2.3 Abstraction with Timing Information 

So far we have used only the sign of the derivative of 3 to determine the switching 
conditions of the hybrid system. However, the time for the system to move from 
one box to its adjacent boxes is omitted, in other words, in the approximating 
system this time can be anything between 0 and +oo. To obtain a more precise 
abstraction, we can include information about the time the system can stay in a 
location, which we call staying time. In order to see how the information about 
staying time can improve the abstraction precision, we notice that in the original 
system a point (x,z) satisfying the condition for the transition from B 1 to B 3 , 
such as g(x,z l ) > 0, does not necessarily lie on the switching surface d(B 1 ,B 3 ) 
and thus from there the system can either continue with the dynamics of B 1 or 
switch to the dynamics of B 3 . Consequently, for soundness, the transitions in 
the approximating hybrid automaton in Figures 1 and 2 are not urgent, and it 
is possible to stay at the same location indefinitely (i.e. no staying condition is 
imposed). It is clear that there may be boxes in which the original system can 
stay for only a finite time; hence, adding constraints on staying time at each 
location allows to reduce approximation error. 

For a general nonlinear dynamics of z, it is not easy to estimate the staying 
time. A method to approximate the smallest time a linear systems with constant 
input can stay inside a convex polyhedron is developed in [15] and extended to 
uncertain linear systems in [11]. However, its generalization to nonlinear systems 
requires solving a nonlinear optimization problem. Here, we propose a simple 
method to exploit timing information by considering not only the sign but also 
the values of i. More concretely, we additionally discretize the derivative of 2 
into a finite number of disjoint boxes, as it is done for z. Each location of the 
approximating hybrid automaton now corresponds to a box B 1 of 2 and a box 
/3 3 of 2 , and we label it with loc 1} . Then, based on the intervals of the derivative 
of each component Z{, we can estimate the bounds on the staying time and then 
embed this information in the transition guards. 
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To facilitate the discussion, we introduce some definitions and notations. Loca- 
tion loc 1J is called a neighbor of location Zoc pq if either the corresponding boxes 
B 1 and B p or the boxes ft and /3 q are adjacent. Location loft is called the left 
Zi-neighbor of !oc pq if Pi — ij = 1 and the right Zi-neighbor if i,; — p, = 1. The 
left, and right ii-neighbors are defined similarly. 

We first specify the staying conditions of the approximating hybrid system. The 
staying condition of location loft is simply z £ ft, which can be rewritten as 
G 3 ( x)C\ft ^ 0 where G 3 (x) = {g(x, z) \ z £ B 1 }. Since g(x, z) is continuous 3 , only 
transitions between neighbor locations are possible. However, we do not know 
the evolution of the derivative of g, therefore we let the guards of the transitions 
between 2 -neiglrbors be true, meaning that the switchings between i-neiglrbors 
are only restricted by the staying conditions of the locations. To define the guards 
of the transitions between 2 -neiglrbors, we shall use the bounds on z to estimate 
the time ty the system can stay within a location loc 13 . Again, we illustrate 
the idea with the case where 2 is scalar, i.e. m = 1. Consider location loc 13 
which corresponds to intervals [z l ,z l ) of 2 and [ft , ft) of z. We distinglruish the 
following 3 cases: (1) If ft < 0, then a transition from loc' 3 to its right 2 -neiglrbor 
is impossible. The staying time is ty £ [{z l — zf )/\ft\, ( z 1 — zf)/\ft\) if ft yf — oo 

and ty < (z‘ — zf)/\ft\ otherwise; (2) If ft > 0, then a transition from loc 13 to 
its left 2 -neiglrbor is impossible and the bounds on ty are defined similarly; (3) 

If ft > 0 A ft < 0, then the transitions from loc 13 to its both left and right 
2 -neighbors are possible. However, unlike in the two previous cases, the staying 
time ty may range from 0 to + 00 . 



2.4 Application to Multi-affine Systems 

We have presented an abstraction method for nonlinear continuous systems. 
The resulting abstract system is simpler than the original system in terms of 
dimensionality, however it requires the ability to deal with nonlinear differential 
inclusions. In the remainder of the paper we focus on the reachability problem for 
a class of differential inclusions which are uncertain bilinear control systems. The 
study of such systems is motivated by our interest in applying the abstraction 
approach to a large class of biological systems which are modeled as multi-affine 
systems [19,23]. Indeed, by projecting out some variables of a multi-affine system, 
one can obtain an uncertain bilinear system, as illustrated in Section 4 where 
we study a simplified model of a biological system. 

Before proceeding with the reachability problem for bilinear systems, we mention 
that besides the interest of bilinear systems for effective applications of our 
abstraction approach, these systems have received much attention over the past 
decades since they could represent a variety of important physical processes 
in engineering. A number of results related to the control of such systems can 

3 In more general cases, the discontinuities in g can be modeled explicitly by discrete 
transitions with resets. 
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be found in [25]. On the other hand, it should be noted that the problem of 
approximating viability kernels of differential inclusions, which is closely related 
to the reachability problem, was studied in [29]. 



3 Reachability Analysis of Bilinear Control Systems 

In this section, we present a method for solving the reachability problem for 
uncertain bilinear systems with both multiplicative and additive control input: 

i 

x(t ) = u(t)) = Ax(t) + y; Uj(t)Bjx(t) + Cu(t ) (4) 

j= i 

where x £ R n is the state variables; u(-) £ U , the set of admissible inputs 
consisting of piecewise continuous functions u of the form u : IR + —>[/,[/ is a 
bounded convex polyhedron in . The matrices A , Bj and C are of appropriate 
dimension. 

The reachability problem for a system with uncertain input can be formulated 
as an optimal control. The essential idea of our reachability method is to use the 
Pontryagin Maximum Principle to find the inputs allowing to derive a conserva- 
tive approximation of the reachable set. 



3.1 Approximating the Reachable Set Using Optimal Control 

Let <p(t, x , u(-)) denote the trajectory of (4) starting from x under input u(-). For 
a set of initial points Ao C and t > 0, the reachable set at time t is defined 
as: lZ(t,X 0 ) = {y | 3 u(-) £ U 3x £ X 0 : y = <p(t,x,u(-))}. Indeed, we can 
show that 7 Z(t,X 0 ) = {x £ R" | V(t,x) < 0} where V(t,x) is the value func- 
tion: V(t,x) = min !i (.) gM {(i 2 (a:o, Ao) | x = ip(t, xo, u(-))} where d(x o, Ao) is the 
distance from xq to Ao- For nonlinear systems, the exact solution V(t,x) can be 
determined by solving a rather complicated HJB equation [20,21]. Reachability 
methods based on solving the partial differential equations have been developed 
(see e.g. [33,24]). As mentioned earlier, our approach is to use the results from 
optimal control to overapproximate the reachable set. More concretely, the idea 
is to track the evolution of the supporting lryperplanes of the initial set under 
some (optimal) input. This idea has been explored in [35,5] to compute polyhe- 
dral approximations of the reachable set of linear control systems. 

Let H be a lryperplane with the normal vector v that supports the initial set A 0 
at point p. Then, for all points x £ Ao we have 

(v,x) - (v,p) < 0 (5) 



where (•,•) is the scalar product. The following result is obtained by applying 
the Pontryagin Maximum Principle (see [27]). 
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Theorem 2. Let S(t ) be the halfspace defined as S(t) = {x £ R" | p(t,x ) < 0} 
where p(t,x) = ( q(t),x ) — (q(t),x(t)) such that q(t) and x(t ) are solutions of the 
following Hamiltonian system with the maximality condition: 

i 

x = Ax + UjBjX + Cu (6) 

i= i 

dH ^ 

q = — —(x,q,u) where H(q, x, u) = (q, Ax + > u^B^x + Cu) (7) 
ox 

j = i 

i 

u(t) £ argmax{(q(t ), UjBjx(t) + Cu) \ u€ U} (8) 

i=i 

wzift initial conditions: q{ 0) = w, i(0) = p. Then, \/t > 0 : TZ(t,Xo) C S'(f), and 
"H(t) = {a; £ | p{t,x) = 0} is a supporting hyperplane oflZ(t,Xo). 

The proof of the theorem can be found in Appendix. We note that the Hamil- 
tonian H in (7) is affine with respect to u, therefore the input u takes its values 
in the boundary of the polyhedron U. Furthermore, we assume the optimality 
of u, and for a bilinear system this assumption can be effectively verified using 
sufficient optimality conditions in [30]. 

Theorem 2 provides a method to overapproximate the reachable set of (4). In- 
deed, by the theorem, for every face of the initial polyhedron Xq there exists an 
input such that tracking the evolution of the face under this input is sufficient 
to derive a polyhedral overapproximation of the reachable set. However, solving 
the optimal control problem (6)-(8) for a bilinear system under a general class 
of input functions is difficult, we therefore restrict to piecewise constant inputs. 
This allows more tractable solutions and, in addition, the error inherent to the 
restriction can be estimated and controled, as we shall show in the next section. 



3.2 Reachability Algorithm for Uncertain Bilinear Control Systems 

Suppose that the initial polyhedron Xq can be represented as intersection of 
a finite number nf of lralfspaces: X 0 = f]"=i where S 1 ' = {a: | {v v ,x) < 
{v v ,p")}, v v is the outward normal vector and p u is a point on the face of S u . 
Let us recall that for tractability purposes, we shall solve the optimal prob- 
lem (6)-(8) in Theorem 2 for a class of piecewise constant inputs. Given the set 
of admissible inputs U and a time step h > 0, we define a set Uh of piecewise 
constant inputs: 

Uh = (u(-) £ U | u(-) is constant on (t k ,t k+1 ), t k = kh, k > 0}. (9) 

We consider the following bilinear equations with v = 1, . . . , n/ which describe 
the evolution of the normal vectors v 1 ' and supporting points p u of the faces of 
the initial set Xq: 
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r f(t)=-(A+E' = 1 ^WBi)Vw, 

I x v {t) = Ax v {t) + Y! j = 1 u’j{t)BjX v (t) + Cu v (t), / 1Q x 

| u u (t) £ arg max {{q v {t),Y^!j=\ UjBjXi(t) + Cu) \ u £ U}, 

[ g"(0) = v v \ x u (0) = p v (initial condition). 

The superscript T denotes the transpose of a matrix. We denote by V(-) the 
polyhedron constructed from the solution (x^ (■) , q" (■) , v," (■)) of (10) under the 
set U of admissible inputs: P{t) = (X£i I (q v (i)> x ) < (^(t), x"(t))}. If 
the set of admissible inputs is Uh defined in (9), we denote the corresponding 
polyhedron by Vhit). Note that by Theorem 2, !Z(t,X 0 ) C V(t). 

Theorem 3. For all t, > 0 the Haussdorf distance between V(t) and Vhit) sat- 
isfies: dH(V(t),Vh(t)) < Ch 2 , where C is a constant depending only on \U\ and 
the norm of the matrices A, B and C of (f). 

The above theorem shows that the error due to the restriction to piecewise 
constant inputs is quadratic in the discretization time step h. This bound is 
proved by using arguments similar to those in the paper [36] , which investigated 
the problem of second order time-discretization of control systems. The proof 
(together with a formula describing the relation between the constant C, |Z7| and 
the norm of the matrices) is omitted due to space limitation and it can be found 
in [11]. In the remainder of this section, we assume that we are provided with a 
scheme to solve the bilinear system (4) under a fixed piecewise constant input 
u(-) £ Ith, which has the form: 

( ^(t fc+1 ) = a v (x v {t k ),u{t k )), 

) q-(t k+l ) = /3’'(q’'(t k ),u(t k )), (11) 

[ t k = kh , v = 1, . . . , nf. 

where u(t k ) is the value of input u(t) for all t £ \t k ,t k+1 ). The development of 
such a scheme is defered to Section 3.3. 

The procedure for overapproximating the reachable set of uncertain control bi- 
linear systems is summarized in Algorithm 1. Each iteration k produces a poly- 
hedral approximation Vh{t k+1 ) of the reachable set at time t k = kh. First, for 
each halfspace S u {t k ) represented by normal vector q"(t k ) and supporting point 
x v (t k ), the value u v (t k ) of the optimal piecewise constant input for the time in- 
terval [t k ,t k+1 ) is computed using the maximality condition. Note that this can 
be done by solving a linear programming problem. Once u v (t k ) is determined, 
the new normal vector q v (t k+l ) and supporting point x v (t k+1 ) are then com- 
puted using (11). Finally, the polyhedron Vh(t k+l ) is the intersection of the new 
halfspaces. Algorithm 1 produces indeed an overapproximation of the reachable 
set of the system (4) with the set Uh of admissible inputs. Using Theorem 3, we 
can enlarge the sets Vh(t k+1 ) by the error bound to obtain an overapproximation 
of the reachable set of the original system. 
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Algorithm 1 Reachability algorithm for uncertain bilinear systems 
for all v £ {!, ,. . . , n/} x l '(0) = v v , q v (0) = p v 
k~ 0 /* (kmax is the maximal number of iterations) */ 
while k < kmax do 
t k+1 = (k + 1 )h 
for all v £ {1, . . . , n/} do 

u v (t k ) £ arg max {(q l '{t k ), YJ l j=i u jBjX v {t k ) + Cu) \ u £ U} 
x v {t k+k ) = a"(x v (t k ),u v {t k )\, q v (t k+1 ) = 0 l/ (q 1 '(t k ),u 1 '(t k )) 
S"(t k+1 ) = {x | (q l '(t k+1 ),x) < {q l '{t k+1 ),x l '{t k+1 ))} 

end for 

Mt k+1 ) = rcLs i '(t k+1 ), k = k + 1 

end while 



3.3 Approximate Solution of Bilinear Systems with Piecewise 
Constant Input 

In this section, we present a method to solve a bilinear system with piecewise 
constant input, which is used as the scheme (11) in Algorithm 1. 

We assume a fixed input u( •) £ Uh such that Vfc > 0 Vf £ [t k , t k+1 ) u{t) = u k £ U 
where t k = kh. For simplicity, we denote by x k = x{t k ) the reachable state at 
time t k under such input u(-). The problem now is to determine the reachable 
state x k+1 from x k . In all the formulas that follow, the superscript k of a term 
is used to indicate that its value depends on the interval k. 

Since the input remains constant during [kh, (k+l)h), given x k one can compute 
x k+1 using the flow of affine vector field A k x + Cu k where A k = A + ^*- =1 u k Bj, 

that is, x k+1 = e hAk (x k + e~ rAk Cu k dr). However, to do so, the transition 
matrix needs to be evaluated for each time interval since A k depends on u k . We 
present in the following an efficient computation scheme which requires matrix 
exponential computation only once. 

The main idea is to consider the bilinear term in (4) during each time interval 
[kh, (k+l)h) as independent input, in other words, the bilinear system is treated 
as a time invariant linear system with input + Cu(t). For 

brevity, we denote W k = Y^j=i u> j ^:i ■ Then, 

x k+i =e Ah x k + f e Mh-r) w k x ^ dT+ f e Mh-r) Cu k dT ( 12 ) 

Jo Jo 

The second integral has a closed form. As for the first integral, we shall approx- 
imate it by replacing the exact solution x(r) for r £ [0, h) with a polynomial 
7 r fc (r) = p k T 3 +p k r 2 +p k r + Pq where p k £ M" satisfying the following Hermite 
interpolation conditions: 7r fc (0) = x k , -?r fc (0) = x(t k ), n k (h) = x k+1 , n k (h) = 
x(t k+1 ). It is well-known that the coefficients of Hermite interpolating polyno- 
mials are uniquely determined [17]. After some straightforward calculations, the 
coefficient p k can be written in the following form: 





42 



E. Asarin and T. Dang 



p k = ( M iU k + Ni)x k + ( P,u k + Q.fix k+1 + r k , i e {0, 1, 2, 3}. (13) 

Then, developing the first integral with x(r) replaced by i r(r) gives: 

3 

x k+ 1 = e Ah x k + r iW k p k + r 0 Cu k . (14) 

2=0 

where R = e A( ^ h ~ T \ l dr , which can also be written in a closed form. Com- 
bining (13) and (14), we obtain an affine relation between x k and x k+1 of the 
form: 

R k x k+1 = D k x k + d k , (15) 

One can see from (14) and (13) that all the terms dependent of u k (i.e. W k 
and p k ) do not involve matrix exponentials. Therefore, using (15) to compute 
reachable states x k+1 , we only need to compute the matrix exponential e Ah . 

Lemma 3. Let x(-) be a solution of (4) under a fixed input u(-) G Uh and x k be 
the approximate solution obtained by the scheme (15) with the same input u(-) 
such that a;(0) = x° . If the derivative x^Ht) is bounded by A4. then for all k > 0 
\x k — x k \ < A4d 4 /(4!). 

The proof of the lemma uses standard results on the remainder term of Hermite 
interpolating polynomials [17], and it is omitted here. The lemma shows that 
the error of the scheme (15) is of order 0(h 4 ). As shown earlier, the error due to 
the restriction to piecewise constant inputs is quadratic; hence, this additional 
error does not change the order of the method. 

4 Application to a Biological System 

In this section, we illustrate our approach with a multi-affine system, used to 
model the gene transcription control in the Vibrio fischeri bacteria. The results 
are obtained using an experimental implementation of the abstraction method 
and the reachability technique for bilinear systems presented in the previous sec- 
tions. This bioregulatory network problem has been studied in [7]. The following 
brief description of the model is borrowed from [7]. The differential equations 
describing the dynamics of a mode of the system are as follows: 



xi = k 2X2 — k\X\Xs + ui 

X 2 = kiXiX 3 — k 2 x 2 (16) 

X 3 = k 2 x 2 — k\X\X 3 — nx 3 + nu 2 



The state variables aq, x 2 , X 3 represent cellular concentration of different species, 
and the parameters k\, k 2 , n are binding, dissociation and diffusion constants. 
The variables u\ and u 2 are control variables (plasmid and external source of 
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autoinducer). We abstract away the variable x\ by discretizing its range and 
construct a hybrid automaton using the abstraction method of Section 2. The 
influence of Xi in the dynamics of the variables X 2 and x 3 is modeled as an 
uncertain input u x , and the resulting system is a hybrid system with uncertain 
bilinear dynamics. As an example, the dynamics of X 2 and x 3 of the location 
corresponding to the interval P of the values of x\ are: 

{ x 2 = kiu x x 3 - k 2 x 2 

x 3 = k 2 x 2 — k\u x X 3 — nx 3 + nu 2 (17) 

u x £ P 

The reachability analysis results for this system are shown in Figure 3. The left 
figure is the reachable set of the approximate hybrid system when the control in 
unactivated (u = 0). The initial set is the rectangle X 0 where x 2 £ [1.05,1.55] 
and x 3 £ [1.25,1.95]. We can see that the uncontrolled system can exit the 
rectangle R = [1, 2] x [1, 2] via the face x 3 = 1 while the control objective is to 
steer the system through the face X 2 = 2. In [7] the following feedback control 
law for this objective is designed: 112 = 6 and U\ is a multi-affine function of 
the state variables x. The reachability computation result for the system under 
this control law from the same initial set Xq is shown in Figure 3 where on the 
left we can see the reachable set of the location corresponding to the interval 
Xi £ [1.0, 1.5] and on the right the reachable set of the location corresponding 
to x\ £ [1.5, 2], This result shows that indeed the controlled system is driven to 
the face X 2 = 2, as desired. 




Fig. 3. Left: Reachable set of the uncontrolled system, i.e. when u — 0. Middle and 
right: Reachable set of the controlled system: location u x £ [1.0, 1.5] (middle); location 
u x £ [1.5, 2.0] (right), the input u x represents the variable xi that is projected out. 



5 Concluding Remarks 

In this paper, we proposed a framework for abstraction of continuous nonlinear 
systems by means of hybridization. We also developed a reachability algorithm 
for uncertain bilinear systems, necessary for an effective application of the ap- 
proach to multi-affine systems. Experimental results are encouraging and various 
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interesting research directions need to be explored. One important issue is the 
choice of variables to project out which still allows to prove the properties of 
interest. In addition, the assume-guarantee ideas from model-checking could be 
useful in this framework. Indeed, one can assume a reachable set for some vari- 
ables 2 which is used as a set of input values in the computation for the remaining 
variables x. The reachable set of x is in turn used as the set of input values for 
the computation on z, which allows to verify the initial assumption. On the other 
hand, we intend to apply our approach to more problems in biological systems. 
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Appendix 

Proof of Lemma 1. We first recall some basic definitions and notations (see [6, 
13] for more details). We denote by B the open unit ball centered at the ori- 
gin. Let T : R ra — > X be a multifunction where X is the set of all nonempty 
compact subsets of R n . T is called upper semi- continuous (respectively lower 
semi-continuous) at x € R™ if for every e > 0 there exists <5 > 0 such that 
Vx' G x + SB tF(x') C T)x) + sB (respectively IF)x) C T(x') + eB). T is called 
one-sided Lipschitz (OSL) with a constant C if and only if for all x,x' G R", 
/ G IF) x) there exists /' G IF) x') such that )x — x',f — /') < C\x — x '\ 2 . 

We proceed to prove Lemma 1. The condition that IF is one-sided Lipschitz is 
easy to verify. To prove that IF is upper semi-continuous, it suffices to prove that 
it is upper semi-continuous at the switching surfaces. Let x be a point on the 
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switching surface S' 1 -*, hence F(x) = F'j(x). For a point x' £ (x+SB) \S^ (in the 
^-neighborhood of x but not on the switching surface), either F(x') = F'(x') 
or F(x') = F*( x'). Since all F 1 are Lipschitz continuous, there exists e > 0 
F 1 (x r ) C F 1 (x) + eB. Therefore, F{x') C F(x) + eB. Obviously, this also holds 
for a point x' £ (x + SB) D hence, T is upper semi-continuous at x. □ 



Proof of Theorem 1. To prove the theorem we suppose that the multifunction 
T in (3) is one-sided Lipschitz with a constant C and bounded on bounded sets. 
We also assume that T is upper semi-continuous and takes closed convex values. 
By Lemma 2, x(-) is absolutely continuous and satisfies dH(x(t),J r (x(t))) < 
A for all t > 0. We shall prove that there exists a solution x(-) of £ 

F(x(t)); x(0) = xq such that for allt > 0 

\x(t) - x(t) | < |x(0) - x(0)|e £ ‘ + e ct - 1). (18) 

To do so, we consider the differential inclusion 

^ £ 0(i(t)); *(0) = x 0 (19) 

where Q(x) = {i> £ F(x) \ { x(t ) — x, x(t) —v)< C\x{t) — x| 2 + A\x{t) — x|}. We 
shall use a well-known existence theorem for upper semi-continuous differential 
inclusions from [12] to prove that (19) has a solution that satisfies (18). To do 
so, we need to verify that Q{x) is nonempty, convex, closed- valued and satisfies 
the condition of boundedness of the solution set. We first prove that Q{x) is 
nonempty for each x. Let w £ F(x) be such that |x(<) — w\ = dn{x{t),J : {x{t))) < 
A. By the OSL condition, given x, we can choose v £ F{x) such that (x(t) — 
x, w — v) < C\x(t) — x\ 2 . Hence, 

(x(t) — x, x(t) — V ) = (x(t) — x, x(t) — w) + (x(t) — X, w — v) 

< \x{t) — x\A + C\x(t) — x| 2 (20) 

The above implies that v £ G(x). It is not hard to see that G{x) is convex, 
closed- valued and satisfies the condition of boundedness of the solution set since 
T does. From Theorem 5.2 of [12] we conclude that there exists a solution x(-) 
of (19). Denoting <j(t) = | x(t) — x(t)|, a is an absolutely continuous function. 
Furthermore, if &(t) exists, then a(t)&(t) = |^cr 2 (t) = {x(t)—x(t),x(t)—x(t)) < 
Ccr 2 (t) + Aa(t), by (20). From this, it can be shown that &(t) < Ca(t) + A for all 
t > 0, and combining with cr(0) = |x(0) — x(0)| we obtain the inequality (18). □ 

Proof of Theorem 2. The equation (7) can be rewritten as q = x, u) = 

— (A+^^- = i UjBj) T q. For brevity, we use x(t) to denote a trajectory <p(t, Xo, u(-)) 
of (4) starting from some point Xq £ Xq under an arbitrary admissible input «(•)• 
Then, d <g (t h a:( * )) = (q(t),x(t)) + {q(t),x(t)) = {q(t), -{A + Y!j=i UjBj)x{t)) + 
{q, Ax{t) + Ej=i UjBjx{t) + Cu(t)) = ( q(t ), - UjBjx(t) + Y ! j=1 UjBjxit) + 
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Cu(t)). By the maximality condition (8), for all admissible inputs u(-) and all 
t > 0, (q(t), + Cu(t)) < ( q(t ), Ej=i Uj(t) B jx(t) + Cu(t )). 

Therefore, < d(q(ty,x(t)) . j n g^dition, f rom the initial condition and (5), 

(g(0),;r(0)) < (g(0),i(0)). Thus, we have Vt > 0 {q(t),x(t)} < (q(t),x(t)). It 
then follows that every point x reachable by (4) from X$ at time t > 0 satisfies 
p(t, x) < 0. □ 
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Abstract. The observability of deterministic, discrete-time, switched, 
linear systems is considered. Depending on whether or not the modes are 
observed, and on whether the continuous state or the mode sequence is 
to be recovered, several observability concepts are defined, characterized 
through linear algebraic tests, and their decidability assessed. 



1 Introduction 

By switched linear systems (SLS), we refer to discrete-time systems that can be 
modeled as follows: 



•Tfe+l — A(^9 k ^)X k T B(^9 k ^)U k /-. \ 

y k = C(9 k )x k , ^ 

where x k € K ra , u k € R m , and y k € are the states, the inputs and the mea- 
surements, respectively. 9 kl which we refer to as the mode at time fc, assumes its 
values in the set {1, . . . , s}, so that the system parameter matrices A(9 k ), B(9 k ), 
and C{9 k ) switch among s different known matrices. We assume that the mode 
sequence { 9 k } ( k Li , whether known or unknown, is arbitrary and independent of 
the initial state and inputs. In particular, we impose no constraints on the time 
separation between two consecutive switches, and we assume that the switches 
are not triggered by state space based events. 

By observability, we mean the ability to infer the initial state x\, and possibly 
a finite portion of the mode sequence (when unobserved) , from a finite number of 
measurements y i, . . . , y^- While the concept of observability has a simple well- 
known characterization in classical linear systems, it has been associated with 
several notions in the SLS literature. Indeed, the fact that the mode sequence 
may or may not be observed, and, in the latter case, that one may or may not 
wish to recover it along with the state, makes for the need to consider several 
different problems, thus different definitions and characterizations. In this paper, 
our aim is to introduce and define several different concepts of observability in 
SLS’s, to characterize them, and to assess their main properties, among which 
decidability is of special importance. 

While we are concerned with the deterministic case, the first attempts at 
defining observability for SLS were motivated by stochastic optimal control and 
optimal estimation problems [7,8,9,11], and are therefore of limited relevance to 



R. Alur and G.J. Pappas (Eds.): HSCC 2004, LNCS 2993, pp. 48—63, 2004. 
(c) Springer- Verlag Berlin Heidelberg 2004 




Observability of Switched Linear Systems 



49 



the present paper. In parallel with those works, the continuous-time case [5], and 
a special discrete-time class [6] were studied, both assuming observed modes. 
The recent surge of interest in hybrid systems expressed by the control and 
computer science communities has motivated renewed research efforts towards 
the SLS observability problem under unobserved modes. Particularly, the work 
in [3] concerned piece-wise affine systems, in which the future mode depends 
on the current state. In [2], two routes were proposed for recovering the modes: 
the first through a finite automaton formalism (which uses extra sources of 
information, beside the measurements, that we, in this paper, assume to be 
unavailable), and the second through failure detection techniques, which require 
the mode sequence to be slowly- varying. Finally, in [10], a systematic attempt 
at carefully defining observability for SLS under unobserved modes was made, 
from which some results in this paper draw inspiration. In the observed modes 
case, observability was characterized and its decidability proven in [1], 

The outline of this paper is as follows. In Section 2, the autonomous case is 
studied. In Section 3, we discuss the non- autonomous case. 



2 Autonomous Systems 



We restrict our attention, for now, to autonomous systems of the form: 

‘Efc+l — A(9k)xk 
Vk — 0(9k)Xki 



obtained simply by removing the B(9^)uk term from (1). Before going any fur- 
ther, we need a few definitions. Since we are dealing with finite-time problems, 
we abandon the mode sequence notation and we define a path 9 as a finite se- 
quence (or string) of modes 9 = 9\9 2 . . . 9n, where N is the path length denoted 
by |0|. We also define On as the set of all paths of length N. Moreover, we 
denote by the infix of 9 between i and j, i.e. ffyj] = . . .9j, we use 

99' to denote the concatenation of 9 with 9' , and we let $(9) = A(9n) ■ • • A(9\) 
denote the transition matrix of a path 9. By convention, we let = e, the 

empty word, and <P(e) = I. We next define the observability matrix 0(9) of a 
path 9 as: 



( 0 ( 9 ,) 



0(9) = 



\C(9 n )A(9 n _ 1 )---A(9 1 ) / 



Finally, we define: 



Y(9,x) = 0(9)x , 



( 3 ) 

( 4 ) 



and we then get that if x = Xt and 9 = 9i9 2 ...9n hr (2), then Y(9,x) = 
(yf . . . yJ ] ) T . Therefore, throughout the remainder of this section, we will use 
(4) to describe (2) in a more compact way. 
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In this section, we will define and characterize several concepts of observabil- 
ity for autonomous SLS (2). The observed modes case is considered in Section 
2.1, and the unobserved modes case in Section 2.2. 

2.1 The Observed Modes Case 

We assume that both the mode path 9 and Y ( 9 , x) are available, and we wish to 
recover the state x uniquely. More precisely, we study the existence of an integer 
N such that the map (x,9) i-»- (Y(0,x),Q) is injective or, in other words, that 

x yf x' => Y(9,x) ^ Y{6 1 x') (5) 

for every path 9 of length N. Since Y(9,x) = 0(9)x, this clearly amounts to 
requiring that 0(9) be of full column rank (i.e. that p(0(9)) = n, where p(-) 
denotes the rank function, and we then say 9 is observable ), in which case, letting 
0(9)^ be a {l}-inverse of 0(9) (see Appendix A), we get 

x = G(9)^Y(9,x). (6) 

The existence of such an integer N was defined in [1] as pathwise observability: 

Definition 1 (Pathwise Observability (PWO)). The SLS (2) is pathwise 
observable (PWO) if there exists an integer N such that every path 9 of length 
N is observable, i.e. satisfies p(0(9)) = n. We refer to the smallest such integer 
N as the index of PWO. 0 

PWO was moreover shown to be decidable. More precisely, the index of pathwise 
observability was shown to be bounded above by numbers Af(s,n ) (actually, 
A f(s,n,n) in [1]) depending only on s, the number of modes (or pairs) of the 
system, and n. Even though the numbers A f(s, n) were not given explicitly, they 
were shown to be smaller than numbers N(s,n) (again, N(s,n,n) in [1]) that, 
moreover, satisfy (we let 1Z(M) denote the row range space of a matrix M): 

Theorem 1 If 9 is a path of length N(s,n), then there exists a prefix 9° of 9 
(i.e. 9 = 9°9 1 for some 9 1 ) and a path 9' of arbitrary length such that 

1l(O(9 0 9')) 01(0(9°)), (7) 

and thus p(O(9 0 9')) = p(O(9 0 )) < p(0(9)). 0 

Even though this theorem was not explicitly stated in [1], its proof can easily be 
derived from that of [1, Theorem 1], 

2.2 The Unobserved Modes Case: Mode Observability 

In this section, we assume that only the continuous measurements Y ( 9 , x) are 
available, and we investigate the possibility to infer a prefix of the path 9 (i.e. 
0njv'i for some N' < |0|) from the successive measurements Y(9,x) only. But 
first, noting that when x = 0, Y (9, x) = 0 for any path 0, we observe that it is 




Observability of Switched Linear Systems 



51 



impossible to distinguish between paths whenever x = 0. As it turns out, this 
happens in general for all states in a union of subspaces of R". Moreover, this 
issue is closely related to false alarms in failure detection, as pointed out in [2]. 
We therefore have to consider the problem from a looser point of view, which 
leads us to the following definition, in which a.e. x stands for “for almost every 
x” , by which we mean for all x € M ra but a union of proper subspaces, thus for 
all x but a set of Lebesgue measure 0: 

Definition 2 (Mode Observability (MO)). The SLS (2) is MO at N if there 
exists an integer N' such that for all 9 € On+n' and, for a.e. x € M n , 

%,AT] ± 0{ hN] =► Y(9, x) ± Y(9' , x') Vx'e R". (8) 

The index of MO at N is the smallest such N' . 0 

In other words, we require the possibility to recover the first N modes (i.e. 0[i,jv]) 
uniquely whenever N + N' measurements (i.e. Y(9,x)) are available, and for a.e. 
state x. To this end, we need a way to discern between the paths 9 using the 
measurements Y(6,x) they produce through Y(9,x) = 0(9)x. As we are about 
to show, the only way to achieve that without any information other than the 
available measurement Y(9, x) is by taking advantage of the following inclusion, 
immediate from Y(9,x) = O(0)x: 

Y(9,x)£*t(0(9)), (9) 

where $t(M) denotes the column range space of the matrix M. The question is 
then whether 9' ^ 9 => Y(9, x) $1(0(9')), which would provide us with a 
simple procedure for recovering a path from the measurements (using the range 
inclusion test, see Appendix A): 

9 = avg 0 , eeN {Y(9,x)e$i(O(9'))} (10) 

The main issue lies in whether the test (10) has a unique solution. In order to 
analyze this, we introduce the concept of discernibility : 

Definition 3 (Discernibility). A path 9 is discernible from another path 9' of 
the same length if 



p([0(9)0(9')])>p(0(9')), (11) 

where [0(9)0(9')] denotes the horizontal concatenation of 0(9) and 0(9'), and 
where the degree d of discernibility is defined as 

d = p([0(9)0(9')))-p(0(9')). (12) 

We then say that 9 is d-discernible from 9’ . 0 

The following proposition is now in order: 



Proposition 1 Y(9,x) ^ $t(0(9')) for almost any x € R" iff 9 is discernible 
from 9'. 0 
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Proof: We have Y ( 9 , x) £ $t(0(9')) iff Y (9, x) also lies in the the output subspace 
of conflict of 9 and 9', defined as: 

0(6,9') =^(0(9))D^(0(9')). (13) 

We therefore need to show that the dimension of the inverse image of C(9 , 9') 
by 0(9), c(9, 9') = 0(9)~ 1 (C(9,9')), which we refer to as the input subspace of 
conflict of 9 with 9', is smaller than n (which implies that its Lebesgue measure 
is 0) if and only if 9 is discernible from 9'. We have: 

dim (C(9, 9')) = p(0(9)) + p(0(ff)) - p([0(9)0(9’) ]). (14) 

Noting that dim(c(0, 9')) = dim(C(0, 9')) + dimker(0(#)), and then recalling 
that p(0(9)) + dimker(0(#)) = n, we get 

dim(c(0, 9')) =n- p([0(9)0(9’)}) + p(0(9')). (15) 

Therefore, by definition of discernibility, we see that dim(c(f?, 9')) < n if and 
only if 9 is discernible from 9' , in which case we moreover have 

dim(c(0, 9')) = n — d, (16) 

where d is the degree of discernibility defined in (12). □ 

Remarks 1 

— Discernibility does not imply that either path is observable (see Example 1). 

— From (16), the degree of discernibility appears as a measure of separation 
between the two paths 9 and 9' : the larger the degree d, the smaller the di- 
mension of the input subspace of conflict. It is clear that the maximum value 
for d is n, in which case the input subspace of conflict is trivial. 

— Note that, as we have defined it, discernibility is not symmetric. 

— jv-i] = 0 m jv_ip *- e - if the two paths only differ by their last value, then 

their index of discernibility is bounded by p. 0 

The last remark raises the question of whether an upper bound (p, the size of 
each measurement yf) is imposed on the maximum degree of discernibility that 
can be guaranteed for all pairs of paths of a certain length. A related limitation 
was pointed out in [10], where p had to equal n, since the criterion considered for 
“switch detection” was similar to ?r-discernibility. It turns out that this limitation 
can be overcome, provided one can use further measurements in order to discern 
the paths, which leads us to the idea of forward discernibility. 

Definition 4 (Forward Discernibility (FD)). Given an integer d > 0, a 
path 9 is forward d-discernible (d-FD) from another path 9' of the same length 
if there exists an integer Nd such that for any pair of paths A and A' of length 
Nd, 9\ and 9' A' are discernible with degree at least d. The smallest such integer 
Nd is the index of d-FD of 9 from 9' . 0 
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Proposition 2 Y(9X,x) ^ ''R(0{9' X')) for all A, A' £ On’ and for almost any 
x £ R" iff 9 is FD (i.e. 1 -FD) from 9' with an index no larger than N' . 0 

Proof: Clearly, the set {x £ R” | 3A, A' £ On’, Y(9X,x) £ Sft(O(0 , A'))} equals 
U A c(9X,9’X '), which, by Proposition 1 and by virtue of the fact that a finite 

union of null sets is a null set, has measure 0 iff 9 is FD from 9'. □ 

We now turn to showing that forward discernibility is decidable. We first 
establish the following lemma, which indicates that the indexes of d-FD increase 
with d, which is not an obvious fact. 

Lemma 1 Let 9 and 9' be two different paths of length N , and A and A' be any 
paths of length N’ . The degree of discernibility of 9X from 9' A' is greater than 
or equal to the degree of discernibility of 9 from 9'. In other words , the degree of 
discernibility is nondecreasing as the length increases. 0 

Proof: It is easily shown, by elementary linear algebra, that 

p([0{ex)0{ffx')\) - P {[0{9)0{9')\) > P (0{9X)) - p{0{9)). (17) 

In other words, the rank of the concatenation must increase by at least the 
increase of each path. □ 

Theorem 2 (Decidability of Forward Discernibility) FD is decidable for 
any degree, as the index ofd-FD, where d is the maximum degree of FD between 
some pair of paths, is smaller than or equal to N(s 2 ,2n). 0 

Proof: Fix 9 and 9', and let A and A' be such that the degree of discernibility 
of 9 A from 9'X' is minimal over all pairs of paths A and A' of length N(s 2 ,2n). 

First, note that the matrices [0{X)0{X')\ are produced by the following set 
of s 2 pairs: 



( A 0 ] A {j) ) ( C (i) C(J)) (i,J) e (1,... ,s} 2 . (18) 

Therefore, by Theorem 1, there exist A 0 and A' 0 , respective prefixes of A and 
X' of the same length, and two paths p and p! of the same, arbitrary, length, 
such that lZ([O(\ 0 p)O(X r0 p')}) C ^([C^A^C^A' 0 )]), which, by [1, Lemma 4] 
and upon some manipulation, implies that 

H([O(6X 0 p)O(9'X' 0 p')]) c TZ([O(9X 0 )O(9' A' 0 )]). (19) 

By Lemma 1, Equation (19) implies that the degree of discernibility of 9X°p 
from 9'X'° p! is equal to that of 9 A 0 from 9’ A' 0 , which, again by Lemma 1, is 
smaller than that of 9\ from 9'X' , which completes the proof since p and p! are 
of arbitrary length. □ 
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Before establishing the main result of this section characterizing mode ob- 
servability, we need the following definition: 

Definition 5 (Complete Forward Discernibility (CFD)). Given an inte- 
ger d > 0, a path 6 is completely forward d-discemible (d-CFD) if it is d-FD 
from every other path 9' of the same length. The index of d-CFD of 9 is the 
maximum index of d-FD, over all 9' ^9, of 9 from 9'. 0 

Theorem 3 The SLS (2) is MO at N iff every path of length N is CFD, (i.e. 
1-CFD). Moreover, the index of MO is the largest index of CFD, over all paths 
of length N . 0 

Proof: Clearly, the set {x £ R" | 3 9,9' £ On, 3A, A' £ On Y{ 9\,x) £ 
3?(O(0'A'))} equals Ue B' Ua,a' c{0X,9'X'). Therefore, by Proposition 2, and by 
virtue of the fact that a finite union of null sets is a null set, it has measure 0 iff 
every 9 is CFD with index at most N' . □ 

We now complete our study of mode observability in autonomous systems 
by answering the following two questions: 

— what effect does N have on MO? In other words, is MO at larger N stronger 
or weaker? 

— Is MO decidable? 

The following proposition answers the first question: 

Proposition 3 If a system is MO at N, then it is MO at any M < N. 0 

Proof: Let N' be the index of MO at N. Then for every pair of paths 9, 9' of 
length N + N' with a switch at or before N (i.e. such that 0i y^ 9\ for some 
i < N), 9 must be discernible from 9'. But, since M < N, this implies the same 
whenever a switch occurs at or before M, which implies MO at M with an 
index smaller than or equal to N' + (N — M). □ 

The converse is unfortunately not true, unless the A matrices are all in- 
vertible (a counterexample is given next): 

Proposition 4 If T(l), . . . ,T(s) are all invertible, then MO at 1 implies MO 
at any positive integer N. 0 

Proof: Let 9 and O' be two different paths of length N, and assume that the 
maximum index of FD over all pairs of different modes (i.e. paths of length 1) 
is N' . It suffices to show that 9 is FD from 9' with index at most N'. Let A and 
A' be any two paths of length N' . It is easy to show that 

N 

c(9X , 9' A') = p| 0(0 [ i,<_i ] )- 1 c(<9 [<jJV1 A, 9^ N ] A')- (20) 

i=l 

Since 9 ^ 9 ' , there exists i < n such that 9i 9[. Note that c{0^ t N] A, Of N ,X') = 
c(9 i p,9 , i p') with p = 0[j +1 ,jv]A, p! = df +1N ,X' and \p\ = \p'\ > N' . Therefore, 
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since, by assumption, 9 t is FD from 9\ with index at most N', dim(c(0[ i> jv] A, 0j,- N , 
A')) < n, and since all the A matrices, and thus are invertible, and 

using (20), we finally get dim(c(0A, d'X')) < n, which completes the proof. □ 

Example 1. Consider 



A(1)= (oo) ■ 4(2)= (oi) (21) 

C-(l) = (10) C(2)= (10). 

The paths 1 and 2 (of length 1) are mutually FD with index 1, but the paths 1 • 1 
and 1 • 2 are not, because they are not discernible and A(l) = 0, which prevents 
further measurements from increasing their discernibility. A 

And finally, 

Theorem 4 MO at any index N is decidable. 0 

Proof: Since the number of paths of length N is finite, and since, by Theorem 
2, FD is decidable, it follows that CFD is decidable, and thus that MO is 
decidable as well. □ 

Note that more precise versions of these last 3 results can be obtained, 
provided one extends the notion of degree to MO. 



2.3 The Unobserved Modes Case: State Observability 

In this section, we are concerned with whether the continuous state x only is 
recoverable. From the previous results, we know that if a system is PWO (with 
index N pwo ) and MO at N pwo with index N mo , then one can recover the state 
x uniquely from Y ( 9 , x) for all 9 of length N pwo + N mo and for almost any x. 
But if we do not need 9 (which is the primary reason behind the “a.e.”), is this 
still the best we can do? It turns out that we can do better, in that we can 
sometimes recover all states x uniquely, for all paths 9 of a certain length. For 
now, we define our concept of state observability: 

Definition 6 (State Observability (SO)). The SLS (2) is SO if there exists 
an integer N (the smallest being the index) such that Vx € M™ and \/9 € 0/v, 

x ± x’ => Y(0, x) ± Y(9', x') V 0' € 6>jv (22) 

In other words, a system is SO if any N consecutive measurements Y(9,x) yield 
x uniquely without knowledge of 9 , i.e. if the map (x,9) >—>■ Y(9,x ) is injective 
in its first coordinate. We first establish a sufficient condition. 



Proposition 5 If a system is PWO with index N pwo , and if every path of length 
Npwo Is n-CFD, then it is SO. 0 
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Proof: Let N c fd be the maximum index of n-CFD, and N = N pwo + N c fd- This 
implies that the dimension of the input subspaces of conflict of any two paths 
of length N satisfying 0[ \,n pwo ] 7^ n pwo } ® therefore the only state 
whose measurements Y(9,x) do not yield 9[i,n pwo ] unambiguously. We then 
have two cases: 

— x ^ 0, in which case the range inclusion test yields 9[i,n pwo \> w bicli can 
then be used in x = 0(9[ 1Npwo ])0}Y(9, x). 

— x = 0, in which case Y ( 9 , x) = 0. By pathwise observability, we then know 

that x = 0. □ 

Example 2. As an example, here is a system satisfying the conditions of Propo- 

sition 5: 

^>=(^)- 4 < 2 >=0l) (23) 

c-(l) = (3 0) C(2) = (2 0), 

Here, N pwo = 2 and N c fd = 2, and it is easy to check that the rank of 
[0(9X)0(9' A')] equals 4 for any pair 9,9' of different paths of length 2 and 
any pair A, A' of paths of length 2. A 

It turns out that the conditions given in Proposition 5 are not necessary (we 
will later give a counterexample). In order to study SO further, we introduce the 
concept of joint observability: 

Definition 7 (Joint Observability (JO)). Two different paths 9 and 9' of 

the same length are jointly observable (JO) if they are both observable, and if 
their left inverses agree on C(9,9'), i.e. 1 

0 0(9 ) {1} - 0(0') {1} )^W) = 0, (24) 

or equivalently , 

(0(9) - 0(9'))P c{efil) = 0 and (0(9) - 0(9'))P c{e> fi) = «• (25) 

Note that, as opposed to discernibility, joint observability is symmetric. A direct 
consequence of this definition is: 

Proposition 6 9 and 9' are JO iff for all x,x' G R", 

x ± x' => Y(9, x) ± Y(9 ' , x'). (26) 

We also need to define forward joint observability: 

Definition 8 (Forward Joint Observability (FJO)). Two different observ- 
able paths 9 and 9' of the same length are forward jointly observable (FJO) if 
there exists an integer N such that for all A and A' of length N, 9 A and 9' A' are 
JO. The index of FJO is the smallest such integer. 0 

1 Given a subspace V, we let Pv denote the matrix of a linear projection on V . 
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Before characterizing SO, we next show that FJO is decidable. 

Theorem 5 FJO is decidable, as the index of JO is bounded by N(s 2 ,2n). 0 

Proof: Suppose that 9 and 8' are observable and that there exist A and A' of 
length N(s 2 ,2n ) such that 9 A and 8' A' are not JO. Similarly as in the proof of 
Theorem 2, we can find A 0 and A' 0 , respective prefixes of A and A' of the same 
length, and two paths p and // of the same, arbitrary, length, such that 

TZ([O(8X 0 p)O(8'X'°p')]) C TZ([O(9X 0 )O(9' A' 0 )]). (27) 

Now, since 9 A and 9' A' are not JO, neither can be 0A° and 9' A' 0 , since Y (9 A, x ) = 
Y(8'X',x') implies Y(9X°,x) = Y(8'X'°,x'). 

Moreover, by Lemma 1, Equation (27) implies that the degree of dis- 
cernibility of 9X°p from 9' A' 0 // equals that of 9X° from 9' A'°0, which 
furthermore implies that c(9X°p, 8'X'°p') = c(0A°, 0'A' 0 ), thus that 

(0(0A» -0(^AV))^c(*Ao f ,,e'A'V) e( i uals (O(0A» -(^'AV^ac.s'A'O) 
and cannot equal zero since its submatrix (0(0A°) — 0(6' X'°))P c ^x° ,0' X'°) is not, 
because, as we have just shown, 9X° and 9'X'° are not JO. Therefore, 9 A°/i and 
9' X'° p' are not JO, which completes the proof since p and p! are of arbitrary 
length. □ 

We now characterize SO: 

Theorem 6 The following are equivalent. 

1. The SLS (2) is SO. 

2. The SLS (2) is PWO with index N pwo , and every pair of different paths of 
length N pwo is FJO. 

3. The SLS (2) is PWO, and every minimally observable path (i.e. a path with 

no observable prefix) is FJO with every other observable path of the same 
length. 0 

Proof: 

2 => 1: Let Nfj 0 be be the largest index of FJO over all pairs of paths of length 
N pW o ■ Let us show that the system is SO with index at most N = N pwo + Nfj 0 . 
Fix a path 9 of length N, and suppose that 9' is such that Y(9,x) = Y(9',x'). 
Let 0pp] be the minimally observable prefix of 9. First, if 8f = 9 [ ipp then 
x = x’ by observability of 0ppp since Y(9,x) = Y(8',x') implies F(0pppa;) = 
Y(8'y 1 x\,x') = Y(9[ i,fc],a/). On the other hand, if QL fc , ^ ^[i,fc]> then since k < 
N pW o , it is easy to show that 9'^ fc j and 0pp.] are FJO with index at most N — k. 
Proposition 6 then concludes that x = x' . 

3 => 2: It is easily seen that the only pairs of paths of length N pwo left to 
check for FJO are those sharing the same minimally observable prefix. Let 9 
and 9' be two paths of length N pwo , and let 9f fc , = 9 pp] be their minimally 
observable prefix. Y ( 9 , x) = Y (9' , x'), which implies Y (0pp] , x) = Y (0' x , x') = 
F^pppa/), implies that x = x' by observability of 0[ip]. 9 and 9' are therefore 
JO, thus FJO. 
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1 => 3: Necessity of PWO to SO is obvious. Suppose that a minimally 
observable path is not FJO with another observable path, i.e. that there exist 
A, A' of arbitrary length such that 0A and 9\' are not JO, which, by proposition 
6, implies the existence of x ^ x' such that Y{9 X,x) = Y{9' A', a/), which 
contradicts SO. □ 

The reason we give two characterizations is that their equivalence is not 
obvious, and because the second one is easier to check, since the number of 
minimally observable paths is in general smaller than s Npwo . Moreover, it is, in 
a sense, much tighter, since two paths can be non FJO only if they do not share 
the same minimally observable path. Finally, 

Theorem 7 SO is decidable. 0 

Proof: PWO is decidable. Since FJO is decidable, and since there is a finite num- 
ber of paths of length N pwo , the first characterization of Theorem 6 concludes. □ 

We now give an example of an SO system that does not satisfy the require- 
ments of Proposition 5: 

Example 3. Let 



(oi) A ( 2 ) ( 03 ) (28) 

0(1) = (10) C(2) = (10), 

This system is PWO with index 2, and any paths of length 2 are FJO with index 
1. For instance, letting 9 = 11, 9' = 22, and A = A' = 1, one gets 



0(9\) 



1 °\ 


/ 1 o\ 




11, 0{9'\') = 


1 2 | , 


(29) 


12/ 


V 1 8 / 





hence that 9\ and 9' A' are JO because the first columns of their observability 
matrices, which span 0(0, 0'), are equal. Thus if we measure Y (/z, x) = (a a a ) T , 
then the initial state can only be ( a 0) T , regardless of the path /./,. A 



3 Non-autonomous Systems 

We now return to the general non-autonomous case, and recall our model: 



Xk + 1 — A{0fc)Xk T JJ(0fc)uZ; 
Uk — O {9}3)x^. 



(30) 
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Our aim here is to extend some of the previous analysis to the system in (30). 
We thus first define, for a path 9 of length N : 



( 



r{6) = 



o 

C(0 2 )B(9 1 ) 

C(9 3 )A(9 2 )B(9 1 ) 



\C(9 n )<P(9[2,n])B(9i) 



0 0 \ 

0 0 



0 



which enables us to further define: 



Y(0,x,U) = O(0)x + r(9)U, (31) 

where U is a control vector in R mA \ Again, if x = Xi, 9 = #i • • • 0jv, and 
U = (uf . . . u t n ) t in (30), then Y(9,x, U) = (yf . . . yJ;) T , and we can concen- 
trate on equation (31). Now, by the separation principle for linear time-varying 
systems, it turns out that we do not need to repeat the analysis of the known 
modes case, and patlrwise observability remains necessary and sufficient for state 
observability in finite time. In this section, we will instead begin directly by tak- 
ing a first look at mode observability. 

Given 9 and 9', our objective in the autonomous case has been, roughly 
speaking, to make the intersection C{9, 9') of Sft(0(0)) with 5 1(0(9')) as small 
as possible. Here, Equation (31) suggests that we should rather consider the 
intersection of the affine subspaces 3ft(C?(9)) + T(9)U and \ R(0(9')) + T{9')U 
(V + v, where V is a subspace and v a vector of R", denotes the affine subspace 
{x + v | x € V}), and study what effect U has on it. Recalling the following 
classic theorem, 



Theorem 8 The intersection ofV + v and V' + v' is either empty or equal to 
V (~l V' + w for some w, in which case it has the dimension ofVC \V' . 0 

we realize that, while the T{9)U terms cannot increase the degree of discerni- 
bility, they can achieve something impossible in the non-autonomous case: they 
can render the output affine subspaces of 9 and 0', i.e. 9?(0(0)) + T(9)U and 
Sft(0(0')) + T(9')U, totally disjoint, which motivates the following definition: 

Definition 9 (Strong Mode Observability (SMO)). The SLS (30) is stron- 
gly mode observable (SMO) at N if there exists an integer N' and a vector U 
such that for all x € R n and all 9 £ 0^+^'; 



%,iv] ^ 9{ hN] => Y(0, x , U) ^ Y{9 ', z', U)Vx'£ (32) 

We refer to such a vector U as a discerning control. 0 

Note that the difference lies in the replacement of “a.e. x” with “V x” , which 
is a stronger statement. In order to characterize SMO, we unfortunately need a 
few more definitions: 
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Definition 10 (Controlled-Discernibility (CD)). Two different paths 9 and 
6' of length N are controlled- discernible (CD) if 

(i - p)(r(9) - r(9')) ^0, (33) 

where P is the matrix of any projection on ift([O(0) 0(9')]). 

It can be verified that CD is well-defined, even though P is not unique. However, 
to fix the ideas, we let P(9, 9’) be the matrix of the orthogonal projection on 
5ft([O(0) 0(9')]), throughout the remainder of this section. Furthermore, note 
that CD is also symmetric. We can now establish the following: 

Proposition 7 If 9 and 9' are CD, then there exists a vector U such that 

Mx £ R", Y(9, x, U) i $1(0(9')) + r(9’)U. (34) 

Even though $1(0(9’)) + r(9')U is an affine subspace, we can still use the range 
inclusion test, by testing whether Y (9, x, U)) — r(9')U is in $i(0(9')). The proof 
of Proposition 7 is as follows: 

Proof: Let U satisfy (I — P(9,9'))(r(9) — r(9'))U y^ 0. Then, by el- 
ementary linear algebra, $i(0(9)) + r(9)U and 3f^(C>(6 ,/ ) ) + r(9')U are 
totally disjoint as affine subspaces of R . pN , which completes the proof, since 
Y(9,x,U) €$t(0(9)) + T(9)U. □ 

Finally, we define: 

Definition 11 (Forward Controlled-Discernibility (FCD)). Two different 
paths 9 and 9' of length N are forward controlled-discernible (FCD) if there exists 
an integer N' such that 9 A and 9'X' are controlled discernible for any pair of paths 
X and X' of length N' . The smallest such integer is the index of FCD. 0 

Unfortunately, we do not know whether or not FCD is decidable. This is in part 
due to the fact that, as opposed to 0(9), we know little about the structure of 
r(9). Nevertheless, we can characterize SMO as follows: 

Theorem 9 The SLS (30) is SMO at N iff any two different paths 9 and 9' of 
length N are FCD. 0 

Proof: Suppose the system is SMO at N with index N’. It follows that there 
exists a control vector U such that for all 9,9' £ On, 9 y^ 9', and A, A' £ On ', 
$i(O(9X))+r(0X)U and $i(O(0' X'))+r(9' X')U are totally disjoint, which implies 
that (I—P(9X, 9'X'))(r(9X) — r(9'X'))U ^ 0, thus that (I-P(9X,9'X'))(r(9X)- 
r(9'X')) y^ 0, hence FCD of 9 and 9' with index at most N’. 

Now, let N' be the maximum index of FCD, over all pairs of different paths 
of length N, and let us show that the system is SMO with index at most N’. 
We need to show the existence of a vector U in R m ( JV+JV 1 such that 



(I - P(9X,9'X'))(r(9X) - r(9'X'))U ± 0 



( 35 ) 
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for all 9, O' £ On, 0 ^ O', and A, A' £ On'- Since every pair 9 and 9' of different 
paths of length N is FCD with index at most N', we get 

(/ - P{9\, 9' A')) (r(0A) - r(0'A')) ± 0. (36) 

Therefore, 

K= |J ker((I-P(9\,9 l \ , ))(r(e\)-r(9 , \'))) ^R m< ^ N+N '\ (37) 

0,8',X,X' 

since it is a finite union of proper subspaces of R m ( Ar + Ar ) ; by (36) . Any control 
vector U £ R. m(N+N )\^ w ju WO rk in (35), and is therefore discerning. □ 

It should be noted that the existence of a single discerning control U im- 
plies that “almost” any vector of the same length is discerning, as established 
by (37). Finally, we describe an SMO system in the next example. 

Example 4- Let 

A(l) = A(2) = (J , (7(1) = (7(2) = (10) (38) 

B(l )=(J) s (!)=(i) (39) 

Since the observability pairs (A(1),C(1)) and (A(2),C(2)) are equal, no two 
paths can be discernible, because all paths of the same length share the exact 
same observability matrix. However, this system is SMO at N = 2, with index 
N' = 1. To see this, it suffices to use Theorem 9 and to establish that ( I — 
P(9X , 9' X'))(r(9X) — r(9' A')) ^ 0 for any two different paths 9 and 9 1 of length 
2, and any pair of paths A and A' of length 1. A 

4 Conclusion 

We have characterized several concepts of observability in switched linear sys- 
tems through simple linear algebraic tests, and we have shown their decidability 
in the autonomous case. An assumption underlying all criteria studied was that 
the mode sequences were arbitrary, which is novel in the sense that most (if not 
all) previous work assumed constraints on the mode sequences, usually in the 
form of minimum “dwell times” between switches. 

This paper is intended as an intermediate step towards a better understand- 
ing of the observability of switched systems. Indeed, some results need to be 
refined, some problems still need to be solved, and many extensions are in view. 
To mention a few, the decidability of forward controlled-discernibility (FCD), 
which seems to be a challenging problem, and the characterization and study 
of state observability in the non- autonomous case, still need to be addressed. 
Finally, the investigation of the application of the concept of discernibility to 
asymptotic observer design promises to be fruitful, and we leave it to a future 
endeavor. 
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A Some Generalized Matrix Inversion Theory 

We now present some definitions and results from matrix inversion theory (see, 
e.g., [4]). We first recall that is a (l}-inverse of O if 

00 {1} 0 = O, (40) 

and that the (Moore-Penrose) pseudo-inverse of O is defined as 

OO^O = O, = e> t , O'O = and OO 1 = (OO^y. (41) 

Note that the pseudo-inverse O * of O always satisfies (40), and is therefore a 
{l}-inverse. If furthermore O has full column rank, then any {l}-inverse 
of O is a left inverse of O, in the sense that O^O = I, the identity matrix. We 
next consider the following equation: 



Y = Ox, 



(42) 
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where x £ R" and Y £ R^, and we examine the conditions on Y for (42) to 
have a solution in x, and how to compute that solution. Note that 

3x\Y = Ox <£=> Y £ 5ft(0), (43) 

which is why we refer to the following test as the range inclusion test: 
Proposition 8 (Range Inclusion Test) If O^ is a {1 ^-inverse of O, then 
Y £ 9t(0) & ( 00 {1} - I)Y = 0. (44) 



Proof: 

4= Let x = O^Y. Then Y = Ox, which concludes the proof. 

=> We have Y £ Sft(O) => 3 x s.t. Y = Ox. By definition of a left inverse, 
we have that OO^Ox = Ox, which implies that OO^Y = Y, which 
concludes the proof. □ 

In words, equation (42) has a solution if and only if ( OO ^ — I)Y = 0 holds for 
some {l}-inverse (it then holds for any (l}-inverse). Note that if (42) admits a 
solution, then x = O^Y is a solution to (42) for any {l}-inverse O W of O. 
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Abstract. We present a parameter inference algorithm for autonomous 
stochastic linear hybrid systems, which computes a maximum-likelihood 
model, given only a set of continuous output data of the system. We 
overcome the potentially intractable problem of identifying the sequence 
of discrete modes by using dynamic programming; we then compute the 
maximum-likelihood continuous models using an Expectation Maximiza- 
tion technique. This allows us to find a maximum-likelihood model in 
time that is polynomial in the number of discrete modes as well as in the 
length of the data series. We prove local convergence of the algorithm. We 
also propose a novel initialization technique to derive good initial condi- 
tions for the model parameters. Finally, we demonstrate our algorithm 
on some examples - two simple one-dimensional examples with simu- 
lated data, and an application to real flight test data from a dual-vehicle 
demonstration of the Stanford DragonFly Unmanned Aerial Vehicles. 



1 Introduction 

The modeling of systems as stochastic hybrid systems has applications in fields 
such as target-tracking, the statistical analysis of time-series data, and systems 
biology. These systems frequently exhibit behavior that is a combination of dis- 
crete switches and continuous evolution; in addition, the data available in these 
applications is usually corrupted by noise. Most target-tracking algorithms for 
maneuvering targets, as well as estimators for hybrid systems, depend on the 
prior knowledge of a good model for the plant dynamics and noise character- 
istics, as well as knowledge of the transition probabilities between the discrete 
modes [1,2]. In this paper, we formulate an algorithm that finds the maximum- 
likelihood values of parameters for both the continuous dynamics in each mode, 
and the transition probabilities between the modes. We draw broadly on sev- 
eral techniques from data association and target tracking [3], motion-capture 
and synthesis methods in computer graphics [4,5] and statistical time series 
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under Software Enabled Control (AFRL contract F33615-99-C-3014) and by an NSF 
Career Award. H. Balakrishnan is supported by a Stanford Graduate Fellowship. 
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analysis [6]. While these techniques are related to classical methods of system 
identification for continuous- and discrete-time systems [7,8], we use them to 
develop a method to identify the parameters of a hybrid system. Given only 
the continuous- or discrete-time output of the system, our algorithm iteratively 
computes the maximum-likelihood parameters for the discrete and continuous 
models, and converges to a local maximum-likelihood autonomous stochastic lin- 
ear hybrid system model. We propose methods to derive good initial conditions, 
so that the local maximum converged to is a suitable model for tracking the 
future behavior of the system. 



2 Model Structure 



We consider a class of hybrid systems, with linear stochastic dynamics in each 
discrete mode. An autonomous discrete-time stochastic linear hybrid system [9] 
is defined to be: 



f x(k + 1) = Aix(k) + Wi(k) 
{ y(k) =Cx(k) + Vi(k) 



k£N 



(1) 



where x £ R" and y £ R p are the continuous state and output variables re- 
spectively. The index i £ {1,2,- •• ,N} represents the discrete state, where N 
is the (unknown, but finite) number of discrete modes in the model. The sys- 
tem matrices are A t £ R nxn for i £ {1,2, ••• ,N} (assumed unknown), and 
C £ R pxn is the measurement matrix ( n and C can be determined using a Sin- 
gular Value Decomposition (SVD) on the output data, and are therefore assumed 
to be known). We denote the covariance of the initial state x(k 0 ) as 7 r 0 € R w , 
and assume that the process noise tUj(fc) and the measurement noise Vi(k ) are 
uncorrelated, zero-mean white Gaussian sequences with the unknown covari- 
ance matrices E [wi(k)wi(k)'] = Qi £ R nxn and E [vi(k)vi(k)'] = Ri £ R pxp 
respectively, where E[-] and (•)' denote expectation and matrix transpose. It is 
assumed that Wi(k) and Vi(k ) are both uncorrelated with the initial state, i.e., 
E[x(fc 0 ) , u:i(/c) / ] = E[a;(fco)^’i(fc) , ] = 0. Y^, T £ M pxT is used to represent the given 
data series, that is, the sequence of p— vectors [y(l),y(2), ■ ■ ■ , y(T)]. We use simi- 
lar notation for other sequences of vectors; for example, the sequence of state vec- 
tors [x(l), x(2), ■ ■ ■ x(t,)) is denoted by Xi- t € M nxt . We denote the set of param- 
eters which defines the dynamics for each discrete mode i by = {A i: Qi, Ri} 
and the entire continuous model by 0 = {0\, 02, • • • , 9n}- 

Given only Y\ .7-, we would like to find a model that maximizes the likelihood ( C 
that the data was generated by this model. To do this, we segment [1,T] into 
the best Ns (> N) segments, such that each segment corresponds to a single 
discrete mode, allowing for modes to be repeated in the data sequence. We label 
segment k as Ik, representing its discrete mode, (Ik £ {1 • • • N}, k = 1 • • • Ns). 
We impose the condition that the system stays in a mode for a given minimum 
dwell time, Tj: this constraint reflects the observation that physical systems do 
not exhibit infinitely fast switching; additionally, we have shown in [10] that 
a minimum dwell time is necessary to estimate the state of a hybrid system, 




66 



H. Balakrishnan et al. 



once the model is given. We also compute the switching times between modes, 
denoted s such that segment k spans the time interval [sfc,Sfc + i — 1], and in 
this segment, the system is in mode Ik- The minimum dwell time constraint 
can be expressed as Sfc+i — Sk > Td- We denote the switching time sequence by 
S = {si, S 2 , • • • sn s }, and the labeling sequence by L = {Zi , I 2 , ■ ■ ■ In s } ( see Fig- 
ure 1). We assume in this paper that the discrete transitions are independent of 
the continuous state of the system, the relaxation of this assumption will be the 
subject of future work. We also assume that the discrete transitions are Marko- 
vian, and we define the Markovian switching matrix M, whose elements are 
M,j = prob(7fc = j\lk-i = *)■ M gives the probability of transition to any mode 
at the switching time. The Markovian assumption is reasonable, since systems 
frequently exhibit probabilistic patterns in their switching behavior - for exam- 
ple, a civilian aircraft is more likely to transition from a turn maneuver mode to 
a straight mode than to another maneuver mode. The optimal segmentation of 
the data sequence into Ns segments can also provide us with the maximum like- 
lihood value of M, so we represent our discrete model by D = {S, L, M}. We use 
1?(.) to denote the parameter space of (•), for example, fi y represents the param- 
eter space of y, and is equal to R p . Finally, we define the function S(statement) 
to be 1, if statement is true; 0, otherwise. In summary, given a data series (of 
length T) denoted by Y\ : t, and knowing the measurement matrix ( C ), we would 
like to find the system parameters, continuous ({A, Qi , I?,}, i = 1 • • • N) as well 
as discrete ({S, L, M}), such that the resultant model maximizes the likelihood 
that the data series was generated by the model. 

For a stochastic linear hybrid system, there are an infinite number of continu- 
ous models that can realize a given output sequence. Since the main problem 
of interest to us is the design of hybrid estimators, we restrict ourselves to the 
class of systems which satisfy the following conditions: (1) We consider the in- 
novations form of the model [7,11], given by x(k + 1) = Aix(k) + Ki(k)ei(k)\ 
y(k ) = Cix(k) + efik), where the innovations e* (k) = y(k) — C'jx(fc), Kfik) is the 
Kalman filter gain and x is the state estimate, and we assume Ai — KiCi is sta- 
ble [11]. (2) We require that the stochastic linear hybrid system be identifiable. 
This condition, while not overly restrictive, roughly ensures that the models are 
distinct enough so as to be distinguished from each other using their continuous 
outputs. This would also give us a class of transformations for the continuous 
model which would explain the given output. The conditions for identifiability 
for a deterministic linear hybrid system have been derived in [12], these can be 
extended quite naturally to stochastic linear hybrid systems. 



3 Parameter Inference Algorithm for Stochastic Linear 
Hybrid Systems 

In this section, we propose an algorithm for hybrid system model inference, assess 
its complexity, and prove its convergence to a local optimum. The structure of 
the algorithm (for one iteration) is given in Figure 2. 
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Fig. 1. Example of a three mode hybrid model, showing model parameters. 



Discrete and Continuous Models from iteration k-1 



rP-v 



Q(k-l) 




Fig. 2. Structure of the Parameter Inference Algorithm for Stochastic Linear Hybrid 
Systems. 



3.1 Expectation Maximization 

For the sake of clarity, we first briefly describe the traditional Expectation- 
Maximization (EM) algorithm [13]. Given general observed data Y £ Gy, we 
could postulate probability density functions (pdfs) g(Y\ip), which depend on 
parameters ip. The Expectation-Maximization algorithm is a method of finding 
a valuation of ip that maximizes g(Y\ip) given the observed Y . Suppose we knew 
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that the observations were in fact incomplete data, for which the underlying 
complete data Z € flz had corresponding pdfs f(Z\ijj). (For example, in the 
case of a stochastic linear system, the observations form the incomplete data, 
while the observations together with the state variables form the complete data) . 
Denoting the vector of parameters by E = (ipi, - ■ ■ ipd), if the mapping from 
flz H > were many to one, we could write 

g{Y\V) = [ f(Z\E)dZ (2) 

To find the model that is most likely to have generated the data, the EM al- 
gorithm iteratively computes log C(E) = log f(Z\E), updating the parameter fit 
(E) at every step. We compute log C(E) instead of C(E), because for many ex- 
ponential families, including systems which are Gaussian, the logarithm of the 
likelihood is easier to compute than the likelihood itself; and since the logarithm 
is a monotonic function, maximizing the log-likelihood is equivalent to maxi- 
mizing the likelihood. However the complete data Z is not available in practice 
(we only have the data we actually measure), and so we use the expectation of 
the log-likelihood, i.e., E^{log C{E) |Y}, derived using the current fit for E. For 
example, in the case of stochastic linear systems, these expectations can be com- 
puted using Kalman recursions. Therefore it is the expectation of the likelihood 
that is maximized; hence the name Expectation- Maximization. 

Algorithm 1 (Classical (G)EM Algorithm [13]) 

Let E^°> be some initial value of E . 

Repeat: 

1. E- Step: Compute 

Q(E\E^) =E* (fc) {log C(E)\Y} = Ey W {log f(Z\E)\Y}, (3) 

where E^(-) denotes the expectation derived using the parameter set E. 

2. M- Step: Choose E^ k+1 ^ to be any value of E £ fly which maximizes 
Q(E\E W), i.e., 

Q(i^(fe+i)|^(fc)) > Q( E\E^), for all EG fly (denoted EM); (4) 
or any value that satisfies the less restrictive condition: 

g(^(fc+i)|^(*)) > g(^(*)|^(*)) (denoted Generalized EM (GEM)). (5) 

until C(E^ k+1 l) — C{E^ k l) converges to e, where e £ R is arbitrarily small. 

The EM algorithm is a special case of the GEM algorithm (5). The EM algorithm 
maximizes the conditional expectation at every iteration, while the GEM only 
ensures its increase. Under fairly general conditions, such as the boundedness 
of the likelihood functions, several properties of the EM and GEM algorithms, 
including convergence, have been proved in [13,14]. In particular, it can be shown 
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that a sequence of likelihoods obtained using the GEM algorithm will converge 
to a local optimum [14]. We will use this property to prove the convergence of 
the algorithm that we propose for stochastic linear hybrid systems, which have 
Gaussian and bounded likelihoods. 

The two main issues with the implementation of any version of the Expectation- 
Maximization algorithm are: (1) the pdf f(Z\F) (or its conditional expectation) 
is sometimes difficult to compute; (2) since the convergence is to a local max- 
imum, a good choice of initial conditions is necessary. We develop methods to 
alleviate both of these problems in the application of EM to the inference of 
stochastic linear hybrid systems. 



3.2 Parameter Inference Algorithm for Stochastic Hybrid Systems 

In order to fit into the form of the standard EM algorithm [13], we note that we 
wish to find the maximum likelihood solution to our model parameters, i.e., 

{D, 0} = arg max £{Yi-.t\D, 0) (6) 

C(Y\-_t\D,0) is analogous to the function g(Y\%jj) in (2). As in the case of the 
classical EM algorithm, rather than work with this density function, we define 
the complete data as the combination of the sequence of state variables (X) 
and the observations (F), from system (1), and denote it by Z. In other words, 
z(k) = {x(k),y(k)} £ R” x R p . The algorithm can be written as follows: 

Algorithm 2 (Parameter Inference Algorithm for Stochastic Linear 
Hybrid Systems) 

Assume an initial continuous model 6b°) and an initial discrete model ZF°b 
Iterate the following until the convergence of the likelihood to a local maximum: 

1 . Step 1 : Find the globally optimal segmentation points (S) and their respec- 
tive labels (L), assuming the model parameters of the current iteration (k). 
Then, update the switching probability matrix M, using: 

Mij = Ylki 2 &{lk - 1 = i)8(lk = j), normalized such that ^ij = 1- 

This gives us the maximum-likelihood discrete model, D^ k+1 h 

2. Step 2: Fit new maximum-likelihood models into the segmented time se- 
quences; i.e., for the computed {S', L} fit the best 0^ k+1 \ 

Drawing an analogy to the classical EM algorithm (Algorithm 1): F = {D,0}. 
Then, we would like to maximize the expectation 

Q(F\FW) = E^ m {logC(F)\Y} = Y DW ^ k) {logf(Z\D,0)\Y}. 



We may rewrite Algorithm 2 in terms of the different conditional likelihoods max- 
imized in each step. The algorithm begins with an initialization of D and 6h°b 
Then, we iterate the following until the convergence of E D (k) qw {log f(Z\F)\Y} 
to a local maximum: 
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In Step 1, we maximize E z yfc) e (fc){fog/(iI|-D,6 )( ^)|E}; this allows us to com- 
pute the conditional expectation E D (fc+i) i6l (fc){log f(Z\0, lM fe+1 ))|y}. 

In Step 2, we maximize E D (fc+i) {log f(Z\0, D^ k+1 ^)\Y}, and compute the 
conditional expectation F, D ( k +i) ^(k+i) {log f(Z\D, 0 < ' k+1 ^)\Y}. 

Algorithm 2 no longer falls into the classical EM framework, since we do not 
compute the expectation in one step and maximize it in the next. In Algorithm 2, 
in Step 1, we assume a continuous model and compute the best discrete model for 
this continuous model. In Step 2, we assume this discrete model, and compute the 
maximum- likelihood continuous model. While both steps correspond to either 
an E-step or an M-step, the expectation computed in one step is maximized in 
the next, and vice versa. Since the likelihood function is changed at every step, 
it is not clear that the convergence properties of classical EM or GEM hold in 
this case, and thus the convergence of Algorithm 2 must be analyzed. 

Theorem 1. Algorithm 2 iteratively generates a sequence of models, whose like- 
lihoods satisfy the model-likelihood sequence conditions (5) of the Generalized 
Expectation- Maximization Algorithm [13,1 f]. Therefore the algorithm is guaran- 
teed to converge to a local maximum. 

Proof. Let us consider the likelihood sequence generated. For iteration k + 1, we 
would like to prove that the model parameters {D^ k+1 \ 6h fe+ b} satisfy (5), i.e., 

Q(£)( fe+1 ), 6>( fe+1) |D (fe ), 0 (fe )) > Q(D (fe) ,0 (fc) |D (fe) ,0 (fc) ). (7) 

We can rewrite Algorithm 2 as follows: 

Step 1: Maximize F D ( k ) 9 (k) {log f(z\D, compute lA fe+1 ) such that 

Q(D( fe+1 ),6>( fc )|D( fc ),6i( fc )) > Q(A6> (fe) |.D (fc ),6> (fe )), for all D e fl D . (8) 
Step 2: Maximize F D ^ k +i) 0 (k){log f(z\0 , D ( - k+1 ' > )\y}: compute 0( k+1 ' ) such that 
Q{D (k+1) ,e {k+1) \D {k \G^ k) ) > Q(D^ k+1 \0\D {k \e ik) ), for all 0 G fl e (9) 
Combining (8) and (9), we get 

S(L> (fe+1) ,6> (fc+1) |D (fe) ,0 (fe) ) > Q(D (fc+1) ,0|D (fe) ,0 (fe) ), for all 0 € fl e 

> Q(D {k+1 \0 (k) \D (k) ,0 (k) ) 

> Q(D {k \0 w \D^ k \0 w ), (10) 

thus proving (7). Since in the case of stochastic models with Gaussian noise, the 
likelihood functions are bounded, the sequence will converge to either a station- 
ary point (saddle surface) or a local maximum. Such a convergence to a saddle 
surface can only occur in the continuous step (Step 2). However, in the case 
of linear Gaussian systems, we can show that the Hessian is always negative 
definite, ruling out convergence to a stationary point. Therefore the sequence 
converges to a local maximum, proving the convergence of the parameter infer- 
ence algorithm for hybrid systems. (Note that the convergence criterion is the 
convergence of the likelihood, and not the model parameters [14]. However, we 
prevent oscillations between different discrete models with identical likelihoods 
by updating the discrete model only when the likelihood increases.) ■ 
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4 Implementation of Parameter Inference Algorithm for 
Hybrid Systems 

Algorithm 2 is of little use if we cannot efficiently compute the likelihoods. In 
this section, we describe the actual likelihood functions chosen, as well as the 
procedures used to maximize them. The problem of maximizing the likelihood in 
Step 1 is potentially intractable, since we need to find the maximum likelihood 
hypothesis from 0(N T ) potential segmentations. However, it is fortunately pos- 
sible in this case to compute the solution, by formulating a dynamic program of 
polynomial complexity [4]. In Step 2, we need to iteratively compute the model 
parameters to maximize the likelihood of the dynamics in each of the modes. 
We use an algorithm which iteratively maximizes this likelihood using a form of 
the EM algorithm [15]. 



4.1 Step 1: Maximizing the Likelihood of the Discrete Model 

Suppose we have initial values of <9 and D (we leave the discussion of details of 
the initialization for Section 4.3). We want to find globally optimal segmentation 
points (S) and labels (L). To achieve this purpose, we employ the following 

dynamic programming algorithm [4]: 

Let us define a reward function £max n (t) as the maximum value of likelihood 
that can be derived from dividing Y\ :t into n segments. £max n (t) is achieved by 
the optimal segmentation of Y\ :t into n parts. Let us define LastMode n {t) and 
LastStart n (t ) to be the label (mode) and start time of the final segment in this 
optimal segmentation. 

Clearly, the optimal segmentation is the one that maximizes the likelihood of 
the entire data set, i.e., has likelihood max K t , £max n (T), and Ns is the 
number of segments in this “optimal” segmentation. Then, for Td < t < T, 

Cmax\{t) = max £(Y 1 . T \6 i ), and LastMode\{t) = argmax£(yi :T |0 i ) (11) 
l<i<N i 

Also, while 1 < n < LyjJ, and nTd <t<T, 

£max n {t) = max [£max n -i(b — l)Mii£(Yb :t Wi)\ 

l<i<N 

(n—\)Td<b<t—Td 



[LastM ode n (t) , LastStart n (t)\ = argmax [£max n -\{b — l)Mn£(Yt, :t \0i)\ (12) 

i,b 

where l = LastM ode n -\{b — 1). In other words, the maximum likelihood de- 
rived from segmenting a sequence into n parts is the maximum product of the 
likelihood resulting from segmenting a smaller sequence into n — 1 parts in the 
best manner, multiplied by the likelihood of the new sequence, multiplied by the 
probability of the corresponding mode switch. A schematic representation of this 
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dynamic programming algorithm is given in Figure 3. Therefore, the required 
optimal solution to the segmentation is 

Cmax(T) = max Cmax n (T) (13) 

and the optimal number of segments is Ns = argma x n Cmax(T). We can also 
find the optimal segmentation ({5, L}) corresponding to this solution. We note 
that 



si = 1; Sjv s +i = T + 1; s n = LastStart n (s n +i — 1), for Ns > n > 1 (14) 
Zjv s = LastModeN s (T)\ l n -i = LastMode n -i(s n — 1), for N$ > n > 1 (15) 

We now also update the switching matrix, using the labels from the optimal 
segmentation 

Mij = 2 = i)S(lk = j), normalized such that M ij = 1- 

This algorithm solves the optimal segmentation problem with a complexity of 
0(NT 3 ), if we know the likelihood function C(Yi :t \0i). In the case of stochastic 
linear systems with Gaussian noise, it is possible to simply express the likelihood 
function in terms of the residuals of a Kalman filter. Given 9i, we can design 
an optimal estimator for the continuous dynamics in the form of a Kalman 
filter. Given the predictions for the continuous state variables E(x(k)\Yi : k~i ) = 
x(k\k — 1), and their covariances P(k\k — 1), we can express the likelihood as : 

I* 1 * 

]ogC(Y 1:t \9i) = --^ lQ g \S ki \ - - YfriWE^nik) (16) 

A— 1 fc= 1 

where r*(fc) = y(k) — Cix{k\k — 1) is the residual at time k, with a covari- 
ance Eki — CiP(k\k — 1 ) C'l + Ri. It can be shown ([13,14]) that maximizing 
C(Y 1:t \G,M) over {G,M} is equivalent to maximizing E[log £(Zj{0, M})\Y], 
where Z is the “complete” data, i.e., the joint likelihood of the observed vari- 
ables (Yi :t , and the state variables Xi :t ). Given the optimal segmentation (5, L), 
the function we would like to maximize is the sum of the conditional likelihood 
functions for the data in each segment. To do this, we update the models in Step 
2 to be the maximum (conditional) likelihood values for the data of each of the 
segments. 



Lmax (t) = max {Lmax n _j(b-l )M U L(Y(b:t) \ 0^) } 
; i,b 



Lmax n _,(b-l ) M u L(Y(b:t)\Q i ) 




l = LastMode ,(b-l) i 
n-1 



Fig. 3. Dynamic Programming Algorithm for Step 1. 
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4.2 Step 2: Finding the Maximum Likelihood Continuous Model 

For each mode, we fit the maximum likelihood model, using a form of the EM 
algorithm proposed in [15,6]. We demonstrate this for one of the segments which 
we (in Step 1) have labeled as mode i. Suppose the data output of this sequence 
is Yb + i : b +n , where b + 1 < b + n < T, n > Tj. Then, the desired log-likelihood 
for this segment can be written as 

11 Tl 

log C(Z b+1:b+n \0i) = -- log |A| - ~{x{b) - n)'E~ 1 (x(b) - n) - - log |Qi| 

, b+n 

-y V, i x ( k ) - Aix(k - l))'Q~ 1 (a;(fc) - +x(fc - 1))(17) 

/c— 6+1 

6+n 

-^logl-Ril - ^ X! (: y{k)-C i x{k))'Rr 1 {y{k)-C i x{k )) 

/c— 6+1 

where, as before, Z b+ \. b + n is the joint (“complete”) data, namely, the observed 
variables Y b+ i ;b+ni and the continuous state variables, X b +i -.b+n'i I 1 and E are the 
mean and covariance of the initial values in that segment of the continuous state 
variable, i.e., x(b). Then, as explained earlier, the maximum likelihood solution 
is the one that maximizes the function 



Q(tt,E,Ai,Qi,Ri) — E\iog£(Zb+l:b+n-l\Oi)\Yb+l:b+n\ (18) 

In computing the conditional expectation in (17), we need to compute the follow- 
ing conditional means and covariances, which are easily obtained using Kalman 
recursions [16]: x fc (fc|s) = E[x(fe)|3'&+i : &+ s ], P b (k\s) = cov(a;(/c)|Yb + i:f>+s), and 
P b {k, k — l|s) = co v(x(k),x(k — l)|Fb+i:b++ We define the following quantities: 

*4 = Efctb+i ( pb ( k ~ !| n ) + x b {k - 1| n)x b (k - l|ra)') 

B = Efctb+i ( pb ( k > k ~ l| n ) + x b (k\n)x b (k - l|n)') 
c = Efctb+i {P b ik\n) + x b {k\n)x b {k\n)') 

As shown in [15], taking conditional expectations on (17), we get: 

Q{n,E,A i ,Q i ,R i ) = -\\og \E\ - |tr {A -1 (P b (b\n) + (x(b) - n)(x(b) - /+)} 
-f log |Qi| - |tr {Qr\C - BA\ - A0 + AiAAQ} 

-f logji?i| 

_ 1 . r D -1 y-b+n [(y(k) - CiX b {k\n)){y{k) - CiX b {k\n))' 

2 1 * Z^fc=6+i + cr < P 6 (jfc|n)C'J]} 



where tr denotes the trace of a matrix. It can be shown, as in ([15]), that at 
iteration r, the choice of parameters that maximizes Q(A i ,Q il Ri) is: 
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(19) 

(20) 



Ai(r + 1) =BA~ 1 

Qi(r + 1 ) = -{C-BA^B'} 

n k J 



Ri(r + 1 ) 



1 

n 



b-\-n 

E 



/c— 6+1 



[(y(fc) - CiX b (k\n)){y{k) - CiX b {k\n))' 
+C t P b (k\n)C '] 



(21) 



Therefore Step 2 consists of the following EM algorithm: For each of the models, 
for the corresponding optimal segment (s) from the E-step, we repeat the follow- 
ing until the estimates and the log-likelilrood converge to a local maximum. 

1. Compute the means and covariances of the continuous state estimates using 
a Kalman filter. 

2. Find the new iterates of parameters using (19), (20) and (21). 

The converged values of the (continuous) model parameters (6*) are the maxi- 
mum likelihood model parameter fits to the segmented data sequence. 

As shown in [17,7], for a large number of samples, the difference between the 
estimate and the true model (local maximum it converges to) tends in distribu- 
tion to a zero-mean Gaussian with a covariance given by ■ Therefore, 

the inverse of the Hessian of the likelihood function gives us a measure of the 
uncertainty of our parameter estimates. 



4.3 Initialization 

We have already mentioned that the EM algorithm is only guaranteed to con- 
verge to a local minimum. The optimal solution is a function of how “good” the 
initial guesses of the parameters are. We present here a novel method of initial- 
izing the continuous parameters. Since we know that the system stays in a mode 
for at least a time Td, we fit in the maximum likelihood model for mode 1 by 
applying the method of the Step 2 described earlier, to Yi-.T d - Given a state esti- 
mate and measurement for the current time step, and a model for the dynamics, 
the Kalman filter gives us a prediction for the measurement at the next time 
step, as well as a region of uncertainty around it, in the form of a covariance 
matrix. We call the ellipsoidal region defined by the predicted estimate and the 
covariance, the validation-gate. We iteratively proceed, one time step at a time, 
and compare the measurement at every time step with the validation-gate from 
the Kalman filter prediction from the previous time step. If the measurement 
falls outside the validation gates of all previously identified modes propagated 
from the previous time step, we initiate a new mode into the system. Doing 
this repeatedly, we estimate the number of modes ( N ), the initial segmentation 
( S,L ), and the initial M. If the stochastic hybrid system is identifiable, then 
the Kalman filter validation gates will be distinguishable, and the algorithm will 
find the correct number of modes; if not, it will only give us a possible model 
that would explain the output data. 
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5 Examples 

Examples 1 and 2: We first present two illustrative 1-D examples to explain 
the proposed algorithm. We use simulated data generated from the stochastic 
linear hybrid system x(k + 1) = CLix{k) + v(k)\y{k) = Cix(k) + w(k) [9] 
where i € {1,2,3}, ai = 0.5, <22 = 1.0, 03 = 1.1, and Cj = 1, for i = 1,2,3. 
The noise covariances are 1.0 and 0.1 respectively. The switching times in the 
original system are at 39s., 56s., and 79s. The results of the parameter inference 
algorithm for hybrid systems for this 1-D example are given in Figure 4. The 
algorithm correctly identifies three discrete modes, and correctly detects the 
switching times. It also converges to a model consistent with the actual model. 
This is illustrated by the following extracts from the algorithm output: 



Parameter 


ai 


<12 


a 3 


Qi 


Q 2 


Q 3 


Ri 


R2 


r 3 


Transitions 


True Value 


0.5 


1.0 


1.1 


1.0 


1.0 


1.0 


0.1 


0.1 


0.1 


[39 56 79] 


Estimate 


0.5014 


1.0564 


1.1023 


1.547 


1.2547 


1.436 


0.0939 


0.0455 


0.156 


[39 56 79] 



We realize that the estimates of the dynamical parameters (at) are more critical 
to this method than the noise covariances, which are statistical, and therefore 
whose values will depend on the number of trials of data available. Clearly, there 
is a very good match in the actual and inferred values of the system dynamics and 
switching times. This is further illustrated by another example, this time with 
data generated by a system with parameters ai = 1.3, <22 = 1.5, <23 = 0.6666, 
and Ci = 1, for i = 1 • • • 3. In this case, with the same switching times as before, 
the system parameters are identified correctly. 



System Parameter 


ai 


a 2 


a3 


Switching Times 


True Value 


1.3 


1.5 


0.6666 


[39 56 79] 


Estimated Value 


1.3 


1.5 


0.6666 


[39 56 79] 



The training data and the results of the mode-detection are shown in Figure 5, 
along with the error plots. The error is the difference between the estimate of 
the continuous state from the inferred model, and the actual continuous state. 
Example 3: We apply the algorithm to data obtained from the DragonFly 
Testbed at the Hybrid Systems Laboratory, at Stanford University. These are 
Unmanned Aerial Vehicles (UAVs), and the data we use are the position and 
velocity estimates. The data shown in this paper corresponds to a Dual-vehicle 
Flight Test. Aircraft 1 is referred to as the evader , and Aircraft 2 as the blunderer. 
The purpose of this experiment that produced this data was to demonstrate 
algorithms for provably safe closely-spaced parallel approaches [18]. However, 
we now use this data for a different purpose: given only the x — and y— position 
and velocity measurements, i.e., no discrete or continuous state information, we 
build a model of the system which estimates the different modes of flight and 
the continuous states. While applying Algorithm 2 to Aircraft 1, we identify 
three discrete modes, with four distinct segments in the data. The switching 
times are estimated to be at t = 9,20 and 39 seconds. Comparing this to the 
times when the autopilot actions were initiated, we find that the system mode- 
transition commands were issued at t = 21 and t = 39 seconds. This further 
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Training data set, inferred modes for 1 D example 1 





Fig. 4. Training data, with output of proposed algorithm for hybrid systems, for ID 
Example (1). Left: Data sets, as a function of time. The different markers correspond 
to different identified discrete modes. Center: Mode transitions inferred from data set. 
Right: Estimation error plot. 




Fig. 5. Left: Training data, with output of Algorithm 2, for ID Example (2). The 
different markers correspond to different identified discrete modes. Right: Estimation 
error plot. 

validates the performance of the algorithm we have proposed. We are also able 
to identify the other segments (and their modes) in the data sequence. In both 
aircraft trajectories, the initial segment (Mode 1) corresponds to the localizer- 
capture segment in which the aircraft first moves into autopilot, and attempts 
to capture the localizer trajectory (straight line). Mode 2 corresponds to the 
localizer-tracking segment. Now that the two aircraft are flying parallel to each 
other, at t = 18 seconds, the blunderer (Aircraft 2) transitions to a mode that 
causes a potential conflict with the evader. Detecting this potential conflict, at 
t = 21 seconds, the evader initiates a maneuver command to avoid the blunderer. 
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Application ol Proposed Algorithm to 
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Fig. 6. Parameter Inference Algorithm for Stochastic Linear Hybrid Systems, applied 
to Dual- Vehicle Flight Test data. Top: (right) Training data, inferred modes, estimates, 
(left) Mode transition sequence. Bottom (left) Error plots, (center) Convergence of 
continuous model, (right) Convergence of hybrid model. 

At t = 39 seconds, the manual pilot takes over from the autopilot and flies the 
aircraft back. The trajectories are plotted in Figure 5, along with the norm error 
of the smoothed estimates compared to the true data, the convergence plots, 
and the mode transition sequences. 



6 Conclusions 

We have proposed an algorithm that determines a maximum-likelihood hybrid 
system model, given only the continuous output of the system. With an intelli- 
gent choice of initial conditions, we also compute physically realistic models. We 
have applied these methods to the analysis of flight data from our UAV testbed 
(DragonFly) and obtained models that allow us to track the vehicles. These 
methods will have a wide-range of applications ranging from system identifica- 
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tion for the many modes of an aircraft, to the problems of modeling biological 
systems as hybrid systems for the purpose of analysis. These techniques also 
provide a means of determining the parameters of a system that are required for 
the efficient design of hybrid estimators, as well as algorithms for multiple-target 
tracking and identity management [2]. 

Acknowledgments. The authors would like to thank Rodney Teo for provid- 
ing the data from the CSPA algorithm verification, as well as for invaluable 
discussions on the same. 
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Abstract. The problem of maintaining the crankshaft speed of an au- 
tomotive engine within a given set interval ( idle speed control), is for- 
malized as a constrained control problem using a hybrid model of the 
engine. The control problem is difficult because the system has delays 
and a large number of constraints. The approach for the synthesis of a 
controller for this system is based on the theory developed for affine sys- 
tems on polytopes. A structured control synthesis procedure is applied 
in which constraints for state and input variables are backward propa- 
gated from the controlled output (the crankshaft speed) across successive 
subsystems. 



1 Introduction 

In the automotive industry, increased performance, safety and time-to-market 
pressure require the use of complex control algorithms with guaranteed prop- 
erties. Best practices in this industry are based on extensive experimentation 
and tuning of parameters for the control algorithm and for the engine model. 
This procedure needs a substantial overhaul to eliminate long re-design cycles 
and potential safety problems after the car is introduced in the market. Using 
more accurate models and control algorithms with guaranteed properties reduces 
greatly the need for extensive experimentation and points to potential problems 
early in the design cycle. In this paper, we investigate this strategy for the idle 
control problem of a four-cylinder engine. 

In particular, we address the problem by proposing 

— a hybrid piecewise affine model of a four cylinder in-line engine that repre- 
sents more faithfully the dynamical behavior of the engine than the tradi- 
tional mean-value models ; 

— the theory of control of affine systems on polytopes. 

* This research has been partially supported by the E.C. project Control and Compu- 
tation IST-2001-33520. 
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The hybrid model for the four-cylinder automotive engine was developed 
by PARADES in close cooperation with Magneti Marelli Powertrain (see ([1, 
2]). The model consists of the series connection of three sub-models: the intake 
manifold, the cylinders, and the crankshaft. 

The difficulty of the problem lies in the load variations coming from the in- 
termittent use of devices powered by the engine, such as the air conditioning 
system and the steering wheel servo-mechanism, which may cause engine stalls. 
The reformulation of the hybrid model as a piecewise affine hybrid system al- 
lows to derive an efficient synthesis procedure for a controller that satisfies the 
specifications. 

The paper is organized as follows: the model is formulated as a piecewise- 
affine hybrid system in Section 2. Section 3 contains concepts and results for 
control of affine systems on polytopes. A control law is formulated for maintain- 
ing the crankshaft speed at a set interval even in the presence of disturbances in 
Section 4. 

2 Hybrid Model of the Engine 

In this section, to set the stage for our control synthesis strategy, the hybrid 
model of the engine described in [1] is cast in the framework of piecewise-affine 
hybrid systems. Piecewise-affine hybrid systems 1 have been proposed in [5,6] and 
are based on piecewise-linear systems as introduced by E.D. Sontag, see [8]. 

Definition 1. A (time-invariant continuous-time) piecewise-affine hybrid sys- 
tem (PAHS) consists of an automaton ( Q , U E c< i, f) in combination with a 

\Q\-tuple of affine systems on polytopes parametrized by q £ Q that interact in 
the following way. At a discrete state q £ Q, the continuous state x q evolves 
according to the affine dynamical system, 

x q (t) = A(q)x q (t) + B(q)u(t) + a{q), x q {t 0 ) = x+, (1) 

y(t) = C{q)x q (t) + D{q)u{t) + c{q), (2) 

with x q £ X q and u £ U. The state set X q for all q £ Q, the input set U, and the 

output set Y are assumed to be polyhedral sets. As soon as a discrete input event 

e £ Ei n is applied, or an event generated by the continuous dynamics occurs, 
e £ E cc i, because the continuous state has reached the guard G q (e) C dX q , a 
discrete transition takes place according to the transition map f and the reset 
map is applied to the continuous state at the past discrete state to yield the initial 
condition at the new discrete state: 

if x~_ = limx g -(s) £ G q -(e) or if e £ Ei n occurs, then 

q + = f(q~,x q -,e), 

x++ = A r (q~,e,q + )x~_ + b r (q~,e,q + ). 

1 The class of piecewise-affine hybrid systems has been proven to be equivalent with 
that of mixed logic-dynamical systems in discrete-time and with linear complemen- 
tary systems, see [3]. 
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Fig. 1 . Engine hybrid model for idle speed control. 

At the new discrete state q + , the evolution of the new continuous state x q + is 
described by differential equation (1), with q replaced by q + , and with initial 
value a: + , . 

9 + 

We recall also the notions of invariant, robust invariant, and controlled robust 
invariant sets that will be used in the rest of the paper. 

Definition 2. A subset of the state set of an autonomous dynamic system is 
said to be forward invariant ( or positively invariant ) if for any initial state in 
the subset the state trajectory for all future times remains inside the state set. A 
subset of a dynamic system with input is said to be controlled forward invariant 
if there exists an input function such that the closed-loop system has a forward 
invariant set. 

A subset of the state set of an autonomous dynamic system is said to be a 
robust (forward) invariant subset with respect to a set of disturbance signals if for 
any initial condition in the subset and for any disturbance signal the resulting 
state trajectory remains inside the subset. A subset of a dynamic system with 
input is said to be a controlled robust (forward) invariant subset with respect 
to disturbances if there exists an input function such that the closed-loop system 
has a robust (forward) invariant subset. 

As shown in Figure 1, the engine hybrid model is composed of three interact- 
ing subsystems, namely the intake manifold , the cylinders and the crankshaft , 
plus the spark and throttle valve actuators. In idle speed control, the output of 
interest is the crankshaft speed n, whose evolution depends on the engine torque 
T, the load torque T) acting on the crankshaft, and the state of the clutch (either 
open or closed). The engine torque T is a function of the spark ignition timing 
and the mass m of air-fuel mixture loaded in the cylinder during the intake 
stroke (the mixture is assumed to be stoichiometric): 

T = r](ip)(Gm + T 0 ) (3) 

where ip denotes the spark advance 2 and rj(ip) denotes the spark ignition ef- 
ficiency. The latter is a strictly increasing function defined on the interval of 

2 The spark advance is defined as the difference between the angular position of the 
crankshaft at the end of the compression stroke and its value at ignition time. 
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feasible spark advances [—15,20], with rj(— 15) = 0.6 and ??(20) = 1. The spark 
advance is positive if ignition occurs during the compression stroke and negative 
if it occurs during the expansion stroke. The spark actuator models the delay 
of the spark actuation system. Due to this delay, the engine control has to issue 
the desired value of spark advance (p at the beginning of the compression stroke, 
for each engine cycle. 

The mass m of air-fuel mixture is controlled by the throttle plate position a 
and is subject to the dynamics of the cylinder filling due to the intake manifold. 

The throttle valve actuator describes the synchronization of throttle control 
with the engine cycle. The throttle valve is driven by sequences of commands 
a , with a time separation of 5 msec, and synchronized with dead-center events 3 : 
each sequence starts with the first command triggered by a dead-center event, 
the number of commands in the sequences depend on the time between two con- 
secutive dead-center events, i.e. on the crankshaft speed. Throttle valve com- 
mands a are nonnegative and bounded from above by 20 degrees. 

Due to the necessary synchronization between the engine cycle and throttle 
valve and spark ignition actuators, the cylinders model returns the dead-center 
event signals to both actuators. 

The hybrid model of the engine has 6 discrete states and a 10-dimensional 
continuous state, i.e. 

Q = {5”_, S, S+, «Sf, S L , S+} and x = (a,T,p,me,m.E,p,PN,T,n,9) 

with: a the throttle valve angular position, t a timer introduced to generate 
throttle valve command events, p the intake manifold pressure, me the mass of 
air- fuel mixture in the cylinder currently in compression stroke, mg the mass 
of air-fuel mixture in the cylinder currently in expansion stroke, <p the spark 
advance, pn the spark advance in the next expansion, T the engine torque, n 
the crankshaft speed, and 9 the crankshaft angle. 

The discrete behavior of the hybrid system, namely the transitions between 
the discrete states and the reset maps, is represented in Figure 2. Transitions 
are triggered by 

— input events in Ei n = {on, off}, modeling opening and closing of the clutch, 
respectively; 

— events generated by the continuous dynamics in E c d = { 9 = 180,0 = 
180 — p,9 = —p, t = 5 } and modeling, respectively, the reaching of a dead- 
center, the actuation of a positive spark advance, the actuation of a negative 
spark advance, and the reading of a throttle valve command. 

Variables a, me, mE,p,PN,T evolve as piecewise constant signals (i.e. according 
to the dynamic x = 0) and are updated by the reset maps defined in Figure 2. 
Note that variables mg and Pn are used only to model negative spark advance. 

3 In 4-stroke 4-cylinder in line engines, at any time each cylinder is in a different 
stroke (intake, compression, expansion, exhaust). Stroke transitions, occurring when 
pistons reach either a top or a bottom dead-center, are synchronous. 
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Consider first the discrete states S, S + and SL , corresponding to the clutch 
in open state. In state S, the cylinder in expansion is generating torque, and the 
cylinder in compression has not received the spark command yet. 

If the spark advance p is positive then, when 9 = 180 — p the spark is 
ignited and a transition to S+ takes place. In state S+ the spark command has 
been given for the cylinder in compression, while the cylinder in expansion is 
generating torque. At the next dead-center, i.e. when 9 = 180, a transition from 
S + to S occurs: 

— the crankshaft angle 9 is reset; 

— the desired value of spark advance ip for the next cylinder is received; 

— the value of the engine torque T produced by the cylinder entering the 
expansion stroke is computed (T + = r](ip)(Gmc + To)); 

— the mass of mixture me loaded by the cylinder at the end of the intake 
stroke is determined (raj = Hp+ M 0 ). 

From state S, if a negative spark advance p is applied then, when the next 
dead-center event is generated (i.e. 9 = 180), a transition to S_ takes place: 

— the crankshaft angle 9 is reset; 

— the desired value of spark advance Pn for the next cycle (related to the 
cylinder currently in compression stroke) is received; 

— the mass of mixture me loaded by the cylinder at the end of the intake 
stroke is determined (mj = Hp+ M 0 ); 

— the engine torque T is reset, since spark has not been ignited yet for the 
cylinder currently in expansion stroke. 

In phase S-, the cylinder in expansion is waiting for the spark command and 
the cylinder in compression has not received the spark command yet. No torque 
is generated in this case. When the specified value p of negative spark advance 
is reached, spark ignition occurs with the transition from S— to S and the cor- 
responding value of engine torque is produced (T + = ff—p)(GmE + To)). 

A timer r, assuming values in the interval [0,5], is introduced to model the 
higher frequency of throttle valve commands. In all transitions produced by 
the dead-center event 0 = 180, the timer is reset and the first value of a new 
sequence of commands is received from the input a. Then, the timer r evolves 
with dynamic t = 1000. The self-loop transitions in Figure 2 detect when r 
reaches the value 5. When they are activated, 5 msec have been elapsed from the 
last timer reset and a new throttle valve command is read from the input a. 

The behavior of the system for the discrete states S^,S L ,S+, related to 
clutch closed, is analogous. Transitions between states S-,S, S+ and states 
St,S L ,S+ are due to clutch switching, which are modeled by the input events 
in E in = {on, off}. 

The model is completed by the continuous dynamics of the p, n and 9. The 
dynamic of the manifold pressure p is identical in all discrete states: 

p{t) = a p (p(t) - po) + b p (a(t) - a 0 ) , (4) 

where po , ao is the equilibrium point around which the linear model has been 
obtained. 
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The dynamic of the crankshaft depends on the clutch state: 

_ f a n n (t) +b n [T(t) -T p - Ti(t)} if q£{S-,S,S+} , , 

[ > \ a^ + bkim-Tp-Tiit)] if qe{S±,S L ,S$} W 

9{t) = 6 n(f) (6) 

where T p is a constant term modeling pumping work and friction, Ti(t) £ [0, T M ] 
is a bounded disturbance modeling variable friction and the action of subsystems 
powered by the crankshaft, n is in rpm and 6 is in degrees. 

Finally, since all the continuous state variables are nonnegative and bounded 
from above, then in each discrete state q the continuous input (tp,a), state 
x and output y = n evolve in polyhedral sets U, X q , Y, respectively. How- 
ever, the boundaries of X q that can be reached, causing an internal event 
which trigger a discrete state transition, are only those defined by 9 = 180, 
for q £ {S, S+, S L ,S^}, 9= 180 - ip, for q £ {S, S L }, 9 = -ip, for q £ {5_, S^}, 
and t = 5 in any q £ Q. 



3 Control Synthesis for Affine Systems on Polytopes 

Our approach to synthesis for idle speed control is based on concepts and the- 
orems derived for affine systems on polytopes. An important role is played by 
the so called control-to-facet problem. The problem is to guide the closed-loop 
trajectory of the system from a given point to a particular facet without first 
crossing or touching any other facet of the polytope. Necessary and sufficient 
conditions for the existence of a control law that solve the control-to-facet prob- 
lem exist [5,7]. We report below some extensions of these results to the case of 
discrete-time affine systems on simplices. 

Problem 1. Consider a discrete-time affine system on a simplex, 

x(t + 1) = Ax(t) + Bu(t) + a , x(to) = xq 

Xi = convh({ui, . . . , u n+ i}) 

with A £ B nxrt , B £ ]R" xrn , a £ fi", and a polytope 

X 2 = convli({u)i, . . . , w m }) = {x £ ]R"|nJx < kj, Vj £ (1, 2, . . . , m}} 

with m > n + 1, normal rij £ ]R" pointing out of the polytope, and Wj opposite 
vertex of normal rij £ E n . Determine a control law g : X\ — > U such that the 
closed-loop system, 

x(t + 1) = Ax(t) + Bg{x{t)) + a, x(t 0 ) = x 0 , 
is such that for all t £T, x(t) £ Xi implies that x(t + 1) £ X 2 . 

A special case of the above problem is that with X 2 = X \ , and then the closed- 
loop system leaves the set Xi invariant or Xi is a controlled invariant set. 
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Proposition 1. Assume that the polytope X 2 is a simplex, i.e. m = n + 1. If 
there exist vectors u\, . . . , u n+ i £ U = IR m such that 

nj (Awi + Bm + a) < kj , Mi,j £ { 1 , 2 , . . . ,n + 1 } , ( 7 ) 



(wf V 



then the matrix 



V^n+l 1 



is invertible and the affine control law g{x) = Fx+h, 



with 




( wf 1 
\ wZ+i 1 



-1 



( u\ 

W+1 



£ ]R( n + 1 ) xrn 



(8) 



solves Problem 1. 



Proof. Since X 2 is a simplex, then the matrix with parameters Wj in (8) is 
invertible, see [7] for a proof. Note that, by (8), for all j £ {1,2 ,...,n+ 1}, 
Fvjj + h = Uj . The closed-loop system is, 

x(t + 1) = Ax(t) + Bu(t) + a = (A + BF)x(t) + (a + Bh). 

Because X 2 is a simplex, 



n+1 n+1 

3c G IR+ +1 , ^ Ci = 1, such that, x G => x = ^ CiWi. 

i = 1 i= 1 

n+1 n+1 

u(t) = Fx(t) + h = CifFwi + h) = '^2 CiUi . 

2=1 2=1 

njx(k + 1) = nJ[Ax(t) + Bu(t) + a] 

n+1 

= ^ Ci nJ ( Awi + Bui + a) < Cikj = kj, V? G {1, 2, . . . , n + 1} 
2=1 

=> a;(f + 1) £ X 2 . 

Note that if X 2 = X\ = X, then X is a forward invariant = subset. If X 2 C + 
then + + 1) £ X 2 C X\, hence x(t + 2) £ X 2 C Xi and X 2 is forward invariant. 



Proposition 2. Problem 1 admits a solution for polytope X 2 not a simplex, if 
there exists a feasible solution for F and h to the problem 

nj ( Am + B(Fvi + h) + a) < kj, Vj £ {1,2,..., to}, Vi £ {1, 2, . . . , n + 1}. 

A technique for computing feedback parameters F and h, based on a partitioning 
of the set X\ is presented in [7]. 
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4 Control Synthesis for Idle Speed Control 

In this section, we use a structured synthesis method analogous to back-stepping 
to synthesize a hybrid idle speed control algorithm. In general, the idle speed 
control problem can be formalized as follows: 

Problem 2. Given a set point for the crankshaft speed n° £ (0, oo) and an inter- 
val I n C (0, oo), such that n° belongs to the interior of I n , determine a control 
law such that the closed-loop system satisfies the following control objectives: 

(1) Bounding the crankshaft speed. There exists a non trivial robust invariant 
set X ri C {x £ X\n £ /„}. Thus, if Xq £ X ri then, for any disturbance 
signals, x(t) £ X ri and n(t) £ /„ for all t > 0. 

(2) Asymptotic stability. In the absence of disturbances, linp^oo n{t) = n°. 

(3) Attraction. For any disturbance signals, lim^oo x(t) £ X„. 

(4) Rejection of disturbances. For stepwise disturbances, lim^oo n(t ) = n°. 

We focus on the control objective (1) formulated above, with I n = [750, 850] rpm. 
The problem is difficult because the dynamic system for the crankshaft speed 
is a system with delays, where the duration of the delays depends on the full 
dynamics. Our approach to controller synthesis makes use of the structure of the 
dynamic system. The system consists of a series connection of 

1. the intake manifold system, 

2. the cylinder system, and 

3. the crankshaft system. 

First, we synthesize a control law for the torque, which is the input to the 
crankshaft system. Then, proceeding backward in the system, we produce a 
control law for the cylinder system with as input the gas mixture and the spark 
advance angle such that the control objectives are met for the torque. Finally, 
we devise a control law for the throttle valve, which is the input to the intake 
manifold system. 

Step 1 - Controller synthesis for the crankshaft system. By specification the 
crankshaft speed n is bounded between 750 rpm and 850 rpm. Furthermore, 
for idle speed control it is reasonable to bound the engine torque to the range 
[0, 30] Nm. Then, consider an extended model of the crankshaft with state vari- 
ables (n,0,T), defined on the multi- variable box 

B c = {(n, 9, T)\n £ [750, 850], 9 £ [0,180],T £ [0,30]} 

= [750,850] x [0,180] x [0,30] . 

Let F° and F e be the facets of B c lying on the subspaces 9 = 0 and 9 = 180, 
respectively, i.e. 

F° = [750,850] x {0} x [0,30] and F e = [750,850] x {180} x [0,30] . 

Consider the control-to-facet problem of determining a subset of B c such that 
F e is the unique exit facet. Since 9 monotonically increases in B c , this problem 
can be stated as follows: 
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Fig. 3. Set S c (on the left) and its partition Si U S 2 U S 3 U S 4 (on the right). 

Problem 3. Find a subset S c C [750,850] x [0,30] such that, under any action of 
the disturbances clutch £ {on, off} and 7](t) £ [0 ,T M ], all trajectories starting 
from (n(0), 9(0),T(0)) = ( no,0,7o ) £ F°, with (no, To) £ S c , reach in finite time 
the exit facet F e and ( n(t),9(t),T(t )) £ B c , until F e has been reached. 

Note that the time required to reach the target set F e depends on the evolution 
of the crankshaft angle 9 , which is affected by both the engine torque T(t) and 
the disturbances clutch and T/(f). A polyhedral set S c for which the following 
proposition holds has been computed: 

Proposition 3. If the initial state (no,0o,To) of crankshaft system satisfies 
9 0 = 0 and (no, To) £ S c , with S c as in Figure 3, then F e is the unique exit 
facet from B c . Furthermore, the worst torque disturbance T)(t ) is constant and 
equal either to 0 or T M and the worst clutch signal is clutch constant and open. 

Let t k denote the sequence of times at which dead-center events are produced. 
From the above proposition it follows that 

Corollary 1. If the engine torque is controlled in such a way that 

(n + (t k ),T+(t k ))eS c , Vfc£^° (9) 

then ( n(t),6(t),T(t )) £ B c for all t > to and the control objective (1) of Prob- 
lem 2 is met. 

To ensure that condition (9) is verified at each dead-center, the crankshaft sys- 
tem is represented by a discretized model that expresses the evolution of the 
system at dead-center times. The discretized model is obtained by assuming 
that the engine torque T and the disturbance torque 7) are constant between 
dead-center times. While assuming 7) constant between dead-center times is 
justified by Proposition 3, this assumption is actually not verified for T in the 
case of negative spark advance. In fact, according to the engine hybrid model, 
the expansion stroke starts in the state S- (or St) associated to a zero engine 
torque and proceeds in the state S (or S L ) with the generation of a constant 






